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ABSTRACT 


The purpose of this thesis is to characterize typical levels of vibrational 
noise on the ocean floor to ascertain the vibration’s effect on possible future 
bottom mounted sensors. The data used for this thesis was obtained from 
publicly available recorded information from four ocean bottom seismometers 
(OBS). The OBSs were located in two geographical choke points: the Luzon 
Strait and west of the Strait of Juan de Fuca. These highly trafficked choke points 
were considered to be a good representation of where these experimental 
bottom mounted sensors might be located should they be built. Unix-based 
seismic processing software available from the Incorporated Research 
Institutions for Seismology (IRIS) proved essential to obtaining calibrated data, 
and the methodology used to get the calibrated data is discussed in detail. The 
results showed that one OBS out of the four was highly variable, with decibel 
levels varying widely from day to day. The other OBSs remained fairly consistent. 
In addition, there were no common discrete frequencies between sensors that 
were in the same geographic area. Recommended future research involves the 
study of environmental effects on the OBSs, additional research to correlate the 
results observed in the Luzon Strait, and a look into the electronic noise floors of 
the OBSs used. 
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I. INTRODUCTION 


A. CONCEPT 

The purpose of this thesis is to characterize typical levels of vibrational 
noise on the ocean floor to ascertain the vibration’s effect on possible future 
bottom mounted sensors. Lawrence Livermore National Laboratory (LLNL) is 
interested in possibly using an experimental sensor on the ocean floor. LLNL 
characterized the noise response for various terrestrial versions of this particular 
sensor, and they expressed interest in collaborating in a study of the vibrational 
noise on the ocean floor, focusing on low frequencies during quiescent periods. 
The data used for this thesis was obtained from publicly available recorded 
information from ocean bottom seismometers (OBS). It is beyond the scope of 
this thesis to characterize all low frequency noise on the ocean floor since noise 
can differ by area and the seismic events occurring during the measurement. The 
data was obtained, then, by focusing on two geographical choke points which 
have co-located OBSs. These highly trafficked choke points were considered to 
be a good representation of where these experimental bottom mounted sensors 
might be located should they be built. 

Discussions were held with scientists from LLNL to find out sensor 
information and, specifically, the frequency regimes, amplitudes, and vibration 
axis that affect the current sensors. The LLNL personnel discussed what an 
underwater version of this sensor might look like, but they were more interested 
in the statistics of the vibrational noise of the ocean floor between large 
amplitude events to determine what to account for in future designs. 

B. RELATED WORK 

A thesis by Jonathon P. Scobo entitled “Ocean Observing Systems” was 
used for this research to start obtaining OBS data and the various seismology 
software to process it [1]. Even though most OBS data is publicly available, OBS 
data can be difficult to access and cumbersome to analyze. Scobo’s thesis 
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proved essential to find the desired OBSs, which seismology computer programs 
to use and which to avoid, and how to request specific time periods of data [1]. 

Most other work that was found with regards to seismic vibrations on the 
ocean floor concentrated on larger seismic events as against the noise 
environment between these events. LLNL expressed less interest in large events 
and more interest in the noise environment during quiescent periods between 
large amplitude events. 

Somewhat related research has been conducted which discusses the 
shortcomings of different types of OBSs due to the bottom type on which they are 
mounted. These shortcomings will be discussed, because the data used for this 
thesis could be affected by these shortcomings. 

C. ORGANIZATION 

This thesis will start with background information discussing the 
Incorporated Research Institutions for Seismology (IRIS) website and the Ocean 
Bottom Seismograph Instrument Pool (OBSIP). The background chapter also 
discusses specific OBS sensors used for this thesis and OBS fidelity on different 
bottom types. The next chapter discusses theory (Chapter III) and includes an 
introduction to seismic waves. Chapter IV, Methodology, discusses the process 
to obtain, process, and analyze the OBS data. Next, the results chapter (Chapter 
V) contains a discussion of the low frequency vibrational noise environment in 
the different areas of interest. Chapter VI discusses the overall conclusions as 
well as recommendations on future work, and the appendices contain the 
MATLAB code and data tables. 

D. OVERVIEW 

Four OBSs in two different geographic locations were chosen and 
subsequently their vibrational noise environments were investigated. Figure 1 
from the IRIS website shows a global view of publicly available OBSs. This global 
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view allows the user to zoom in on and select the desired OBS, which then 
opens a link to download the different types of available data. 



Figure 1. The Global View of Publicly Available OBSs. Source: [2]. 

Two geographic choke points were chosen: the Luzon Strait and west of 
the Strait of Juan de Fuca. IRIS publicly available software, specifically a 
program to read information in the form of Standard for the Exchange of 
Earthquake Data (SEED), called “rdseed,” a program used to evaluate response 
files called “EVALRESP,” and Seismic Analysis Code (SAC) were used to 
transfer downloaded data into a usable format and then calibrate the data. The 
analysis conducted on the calibrated data was accomplished with MATLAB. 
Figure 2 and Figure 3 display the available OBSs at the chosen choke points. 
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Figure 2. Available OBSs in the Luzon Strait. Source: [2]. 



Figure 3. Available OBSs West of Strait of Juan de Fuca. Source: [2]. 

The data was obtained, then, by focusing on known geographical choke 

points which have co-located OBSs. These highly trafficked geographic choke 

points were considered to be a good representation of where these experimental 
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bottom mounted sensors would be located for either testing or operations after 
they are built. Initially, more geographic choke points with OBSs were inspected 
in hopes to add additional locations to analyze. One of the OBS locations near a 
choke point south of the Bering Strait was inspected, but the OBSs were short 
period OBSs (discussed later) that were used during an active sound 
transmission test. At this location, there was only a week’s worth of noisy data. 
According to IRIS in Figure 1 there are no OBS locations for the majority of the 
operationally significant choke points. The OBS locations are understandably 
chosen for geological significance (e.g., hydrothermal events in the middle of the 
Atlantic Ocean). 
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II. BACKGROUND 


A. IRIS AND OBSIP 

All of the OBS data obtained for this thesis was from the IRIS and OBSIP 
websites. IRIS is a non-profit organization, created in 1984 by the National 
Science Foundation (NSF). According to their website, IRIS involves over 100 
universities in the United States. Their programs “contribute to scholarly 
research, education, earthquake hazard mitigation, and verification of the 
Comprehensive Nuclear-Test-Ban Treaty” [3]. One of IRIS’s mission statements 
is to promote seismic data availability, and they allow access to data by using a 
simple request format, which will be addressed later in this thesis. 

IRIS oversees and coordinates multiple seismic networks, which includes 
the network used for this thesis, the OBSIP. OBSIP focuses on marine 
seismology, geology, and geodynamics. The OBSIP website [4] also contains 
links describing current OBS instrumentation, links to access OBS data, and 
information for researchers to request OBSIP instrument usage. The OBSIP 
allows their OBS equipment to be available to NSF-sponsored investigators and 
to other government research or educational institutions [4]. 

B. OBS 

Ocean Bottom Seismometers can be broken into two fundamental types: 
long period and short period. According to OBSIP’s website [4], the long period 
OBSs have a 3-component broadband seismometer (one vertical component and 
two horizontal) and a low-frequency differential pressure gauge 
(DPG)/hydrophone. These long period OBSs are primarily designed for passive 
recording and are normally capable of a greater than 12-month deployment. The 
short period OBSs consist of either a vertical seismometer or, like the long period 
OBS, three orthogonal seismometers and a hydrophone [4]. The OBSIP website 
[4] also states that the sampling rates for the short period OBSs are typically 
much higher (200-250Hz), and these instruments are normally deployed from 1 
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to 6 months. The short period OBSs are typically used for active-source 
experiments [4]. 

The OBSs used in this thesis are long period. Because of how OBSs are 
placed on the ocean floor, the orientation of the horizontal channels are usually 
unknown. Methods to orientate the horizontal channels can involve analyzing 
data from the known location of an active source (e.g., air gun) or using P or 
Rayleigh waves from a known earthquake [5]. The orientation of the horizontal 
channels used in this thesis were unknown. 

1. LDEO OBS Sensor MK2 

Figure 4 from the IRIS website displays the highlights from the Lamont- 
Doherty Earth Observatory (LDEO) Standard Ocean Bottom Seismometer. The 
LDEO OBSs used for this thesis are located in the Luzon Strait and were used in 
the Taiwan Integrated Geodynamic Research (TAIGER) project. The TAIGER 
project is a joint USA-Taiwan program that studies the tectonic development of 
Taiwan [6]. The network identifier “YM” designates the OBSs in the TAIGER 
project. The two OBS stations that were used to provide data were YM 38 and 
YM 39, displayed in Figure 5. 
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Ooeaa*BottoRi Sdimology Laboratory 

LDEO Standard Ocean-Bottom Seismometer ' ■> ' i- i = 



The LDEO standard seismometer 
design has been in use for nearly 10 
years. The LDEO OBS lab has built 
and operates 25 standard OBSs as 
part of the NSF OBS Instrumentation 
Pool. The seismometer sensor is an 
L4C 1 Hz geophone, with a low-noise 
amplifier, giving useful response 
down to 100 s, and a differential 
pressure gauge. This instrument has 
been used in both year-long passive- 
source and shorter-term active- 
source experiments. The design 
includes dual redundancy with two 
transponders and two dropweights. 

Variants and add-ons 

Trawl-rcsistant OBS 
Moored DPG/hydrophonc 
Ocean-bottom magnetometer 
DifTusc flowmeter 
Trillium Compact sensor 
Absolute pressure gauge 
Hydrophone 


Specifications 


Max. Depth 

5000 m 

Max. Duration 

400 days @ 125 sps 

Channels 

4,24-bit recording 

Sensors 

L4C 3-component geophones; 
differential pressure gauge 

Response 

100 s - 60 Hz (seismometer) 

0-20 Hz (DPG) 

Leveling system 

Active 360°,motor-driven 

Weight 

750 lb in air 

Footprint 

3'X4' 

Flotation 

9X12'glass spheres 

Sampling 

40-100-125 sps 

Release 

Dual dropweights 

Acoustics 

Two ORE 12 kHz transponders 

Power 

Lithium battery pack, +/- 7.5 V 

Oscillator 

Seascan 10 MHz clock 

Sensor housing 

17’glass sphere 

Burnwires 

LDEO design 

Recovery aids 

Radio, strobe, flag 

Recording 

2 X 32 Gb CompactFlash cards 

Dropweights 

Two steel weights (75 lb in air) 

Datalogger 

LDEO ultra-low power OBS 
datalogger (300 mW @) 125 sps) 


Deployment History 



Figure 4. Highlights of the LDEO OBS, Which Includes a Useful 
Instrument Response Down to a Period of 100 Seconds. 

Source: [7]. 
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Figure 5. Location of YM38 and YM39, the Two OBS Locations in the 

Luzon Strait. Source: [2]. 

The LDEO OBS uses an L4C 1 Hz geophone as the seismometer sensor 
[8], which is displayed in a drawing in Figure 6. As discussed in [9], the vertical 
L4C geophone produces a voltage which is proportional to the difference in 
velocity between the mass, M, and the ground The mass, M, is 

suspended with a spring with stiffness k and a dashpot, b, for damping. The 
voltage is produced by the motion sensing coil with inductance, U, and coil 
resistance, Rc. The L4C 1 Hz geophone has a harmonic frequency of 1 Hz, and 
the resistor, Rs, is an external resistor to shunt the output and critically damp the 
system at that frequency, as explained in the paper by Bowden [9]. 
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Figure 6. Drawing of Vertical L4C 1 Hz Geophone, Which Is Used as the 
Seismometer Sensor of the LDEO OBS. Source: [9]. 

2. SIO ABALONES OBS System 

The Scripps Institution of Oceanography (SIO) Autonomous Broad 
Application Low Obstruction Noise Exempt System (ABALONES) OBS system 
can have a variety of configurations. The particular OBSs used for this thesis, 
located to the west of the Strait of Juan de Fuca as displayed in Figure 7, are 
considered an intermediate-period sensor because of their large frequency 
range. In this configuration they concentrate on long-period vibrations with a 
Nyquist frequency of 25 Hz. Table 1 displays the specification for the SIO 
ABALONES system, and Figure 8 displays a 3D-CAD drawing of the ABALONES 
OBS design. The yellow shell of the ABALONES design is a trawl and current 
resistant enclosure, allowing these OBSs to be set in shallow, continental shelf 
areas to as deep as 6000m [10]. It is equipped with a Nanometrics Trillium 
Compact Seismometer, which has the advantage of a nearly flat response curve 
from 120 seconds to 100 Hz [10], [11 ]. 
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Figure 7. 7D J73A and 7D M02A, the Two OBS Locations West of the 
Strait of Juan de Fuca Used for this Thesis. Source: [2], 


Table 1. SIO ABALONES OBS Specifications. Source: [11 ]. 


Seismometer/ 

Accelerometer 

Intermediate-Period: Nanometrics Trillium-Compact seismometer, which 
has a flat velocity response from 120s to 100 Hz. The seismometer is 
housed within an 8” Inside Diameter (ID) Titanium cylinder and a custom- 
designed active leveling system which can be leveled from ±180° to 
vertical with ~0.3° accuracy. 

Pressure Sensor 

Short-period: Customized HiTech hydrophone with internal preamp. 
Bandwidth (-3 dB) is 50 mHz - 15 kHz. 

Long-period: Differential Pres.sure Gauge (DPG) with response from 10 
mHz - 10 Hz. 

Digitizer 

24-bit A/D with solid-state recording (CompactFlash). The dynamic range 
is ~ 126 dB (3-bit self-noise). The ABALONES design includes the latest 
seismic A/D chip with programmable sample rates from 1 Hz to 4 kHz. 

Clock 

Low-power, digitally-temperature-compensated (DTCXO), precision time 
base. The drift rate is 1:3-5 x 10* (<5 ms/day before correction and <1.5 
s/yr) 

Recording 

Capacity 

Single CompactFlash card slot -> 64 Gb CF cards currently available. 
Capable of recording 1-year deployment on 4-channels @ 100 Hz. 

Data Offload 

USB2 through the end cap without opening the pressure case. Offload 
rates are CF card dependent, currently up to 400 Mb/s. 

Battery Pack 

Batteries are mounted in the main data logger pressure case and optionally 
in an additional battery case. Alkaline batteries power the acoustic release. 

Recording 

duration 

IPA: Lithium battery pack can provide power for 12 months. Alkaline 
packs can provide power for 4+ months. 

SPA: Lithium battery pack can provide power for 18+ months. 

Weight 

850/550 lbs with/without anchor (in air) 

Pressure Case 

5” diameter A1 cylinder - 6 km depth rating 

Release 

Double bum wire operated acoustically 

Dimensions 

29" high X 42" x 42" (w/ bail) 

Power (Total) 

Intermediate: <300 mW (4 chan), plus ~160 mW w/ Trillium-Compact 
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Figure 8. 3D-CAD Drawing of ABALONES Ocean Bottom 
Seismometer Design. Source: [11]. 


C. OBS FIDELITY 

As mentioned by Sutton et al., “The usual reason for placing a 
seismometer anywhere is to measure free-field particle motion in response to 
some seismic source” [12]. The complication presented by soft sediment of the 
ocean floor is that the coupling between the seismometer and the sediment may 
be unknown [12]. Vertical vibration, channel Z in this thesis, is of more concern 
than any horizontal motion, channels 1 and 2. Duennebier and Sutton maintain 
that vertical motion across the water-sediment boundary is reasonably 
represented by an OBS due to particle motion continuity across the interface 
[12], [13]. The fidelity of horizontal motion captured by an OBS is also less 
accurate because of the shear discontinuity between the sediment and the water 
[14]. Problems with OBS horizontal fidelity are compounded by bottom ocean 
currents in addition to the soft sediment. 
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In [14], Duennebier and Sutton point out that horizontal motion is 
discontinuous across the boundary of the ocean floor. This differential motion 
generates torque on the instrument which in turn causes both horizontal and 
vertical motion to be captured by the seismometer channels. In addition, if the 
sensor is not completely level, components of vertical and horizontal motions will 
couple into the wrong channels. Consequently, cross-correlation between the 
horizontal and vertical channels is common on bottom mounted OBSs [14]. The 
horizontal channels are considered to be more susceptible to noise from this 
source. Duennebier and Sutton collected data in a study proving increased OBS 
fidelity when the OBS is buried beneath the sediment on the ocean floor [13]. 

Any sensor placed on the ocean floor will be subject to the same 
vibrations as the OBS. Therefore, these issues relating to the fidelity of the 
seismic measurements may be of concern to LLNL for their potential to affect the 
performance of the proposed sensor as well as for possible limitations in the 
accuracy of the seismic measurements provided in this thesis. If a bottom- 
mounted sensor is used, [14] explains that the effects of unwanted vertical 
motion can be minimized by increasing the footprint of a sensor package relative 
to its height. The tradeoff is a decreased sensitivity to certain types of waves 
when the sensor base approaches too high a fraction of the seismic wavelength, 
particularly in higher frequency shear waves. 
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III. THEORY 


Seismic waves can be divided into two categories: body waves and 
surface waves. The body waves travel through the interior of the earth and are 
the first waves to arrive after a seismic event [15]. The surface waves cause the 
most damage associated with an earthquake, and they travel along the earth’s 
crust [15]. The mathematical theory behind both types of waves is discussed in 
this chapter. Figure 9 is the coordinate system used in the equations. All of the 
equations in this chapter are from Bullen and Bolt [16]. 


X3 


Xl 


Figure 9. Coordinate System Used in Equations 1 through 20. 



A. BODY WAVES 
1. P Wave 

The P wave is the primary or “push wave.” Bullen and Bolt [16] describe 
the P wave in terms of the dilatation, 6. The dilatation, or irrotational disturbance, 
is the fractional change in volume of a point in the medium as the wave passes. It 
is given by the divergence of the vector w, , which gives the displacement of a 

point with respect to its equilibrium position in direction x,. as shown in Equation 

1. Equation 2 is the scalar wave equation for the dilatational disturbance. The P 
wave is transmitted through a substance at speed a , given in Equation 3. 


^ du, 5m, 5m, 

9 = —- + —^ + ^ 

5xj 5^2 5^3 


( 1 ) 
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(2) 


d^e 

de 




a - 


P ) 


( 3 ) 


A is Lame’s first parameter, which is an elastic modulus defined by 
Equation 4. ju is Lame’s second parameter, which is also known as the shear 
modulus, and p is the volume density of the material. The Lame parameters 
together characterize the linear elasticity for homogeneous isotropic media [17]. 
The relationship between Lame’s first parameter and the shear modulus, p , is 
given by [16] 

A = K-^p = p{a^-ip^), (4) 

where K is the bulk modulus of the material. 


Figure 10 illustrates how a P wave travels. It is a longitudinal wave 
consisting of changes in the density of the material. 


T = 0 


T = 1 


T = 2 


T = 3 


Figure 10. Illustration of P Wave. Source: [15]. 
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S Wave 


S waves are a type of shear wave and are described by the vector wave 
equation 


dt^ 


curl(u^) = P^V^curl{u^) 


( 5 ) 


where /3 is the speed of the rotational disturbance [16] and given by 


P = 


\p) 


( 6 ) 


Since the speed depends on the shear modulus, S waves will not 
propagate through fluids which cannot support shear stress. The S wave in the 
far field can be plane polarized [16]. If an S wave is traveling horizontally and the 
particle motion is also horizontal, it is denoted SH. If the particle motion is vertical 
only, the S wave is denoted SV [16]. Figure 11 illustrates an SV wave assuming 
that z is vertical and x and y are horizontal. 


17 



T = 0 

T= 1 

T = 2 

T = 3 


Figure 11. Illustration of an S Wave. Source: [15]. 

B. SURFACE WAVES 

In Figure 12, which is adapted from an illustration in Bullen and Bolt [16], 
M is an elastic, homogeneous half-space below the plane denoted by the Xj and 

^2 axes. M is bordered above the plane by another homogeneous elastic 
half-space, M’. 




Figure 12. Two Isotropic Flomogeneous Elastic Media M’ and M. 

Adapted from [16]. 
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Waves can travel along the boundary between these two half spaces. To 
be trapped along the boundary, the amplitude of these waves must decrease 
exponentially as the distance from the boundary increases. Consider such a 
plane wave traveling in the Xj direction. Since it is a plane wave, the particle 

displacement along any line parallel to Xj is equal. As a result, any partial 
derivative taken with respect to X 2 will be zero [16]. Letting ^ and y/ be functions 


which can be used to express the displacement components and as given 
by 


d(j) dy/ d(j) dy/ 

5 x, 5X3 5X3 5 x, 


( 7 ) 


it can be shown that 


^2 Su, du. 

VV = —-—^ 

5X3 5 x, 


( 8 ) 

( 9 ) 


This suggests that ^ is linked with P waves since they are described by 
the dilatation, y/ is linked with SV waves since it is associated with a differential 
displacement between the direction of propagation and the vertical. Any 
displacement, would be linked with SH waves. 


The wave equations associated with ^ ,y/ , and u,^ are given by 


ay 

df 


= V 


d^y/ 


= V 
= 


( 10 ) 

( 11 ) 

(12) 


where the speeds in M’ are replaced with a' and P'. 
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To find solutions for y/, and [16] sets r and 5 to the following, 

where c is the speed of a set of simple harmonic waves with wavelength Inik . 
Again, for M’, the variables are r' and s'. The variables r, r', s, and 5' 
designated in Equations 13 and 14 are imaginary. 


r = 


f 2 


c 



V 

1 ,5 = 

) 


C 2 


C 



1 1 


2 A 

C 

2 

f 2 \ 

c 

-T-1 

,5’ = 

—T-1 

U" J 


[p'" J 


Solutions in M and M’ are sought in the form: 

^ _ J^^k(-rxy+x^-ct) 
yy — Q^ik^-SX^+X^-Ct) 

U = 

^ ^ik(-r'x^+Xi-ct) 

yy — ^ik{-s'x^+X^-Ct) 
jy _ p\ ^ik{-s'x^+x^-ct) 


(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 


Since both r and 5 are imaginary, the positive or negative square root is 
taken as needed to ensure that all solutions consist of waves that exponentially 
decrease as the distance away from the boundary is increased. Expressions for 
the displacements and as a function of time and position can be found from 
these solutions by applying the relationships which implicitly define ^ and y/ in 
Equation 7. 

The boundary conditions dictate which solutions propagate. For example, 
if M is not homogeneous elastic half space, but is instead a vacuum (or air) the 
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resulting wave is a Rayleigh wave, where there are no SH waves. The particle 
motion is retrograde, and displayed in Figure 13. 


T = 0 


T= 1 


T = 2 


T = 3 


Figure 13. Illustration of Rayleigh Wave. Source: [15]. 

If M is not a homogeneous elastic half space, but is instead an 
incompressible fluid (water) the resulting wave is a Scholte wave [18]. The 
particle motion again is retrograde, and the amplitude decays exponentially with 
distance from the boundary. An illustration of a Scholte wave is displayed in 
Figure 14. 




Figure 14. Illustration of Scholte Wave. Source: [19]. 
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IV. METHODOLOGY 


A. OVERVIEW 

This chapter explains the methods chosen to calibrate the data using IRIS 
software, why these methods were chosen, and the data analysis conducted with 
MATLAB. Multiple routes to calibration were available, and picking the most 
effective way forward proved challenging. Many different software calibrations 
methods are available from IRIS, and a few different methods to get calibrated 
data on MATLAB (manually) were attempted before the steps were decided 
upon. The steps laid out in this chapter should allow follow-on study with 
somewhat less of a learning curve for non-seismologists. 

B. CHOOSING THE OBS 

The global view on the IRIS website, displayed in Figure 1, of the available 
OBSs proved a convenient place to start the OBS selection process. This global 
view allows the user to focus on the desired area and select individual sensors. 
Once selected, a window, which is shown in Figure 15, is opened that displays to 
the user the dates of available data and a link to more information. It becomes 
readily apparent if the selected OBS is a short period vice long period. In most 
cases, the short period OBS dates available will be for days and weeks versus 
multiple months. 
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Figure 15. Convenient Use of Global View for OBS Selection. 

Source: [1]. 

Once the user clicks on the “more information” link, the specific OBS page 
opens and allows the user to view the network and station information with the 
channels available. Figure 16 shows this page. The user can select any of the 
channels and view the response curve, the depth, sensor location in latitude and 
longitude, and the response (RESP) and poles and zeroes (PZ) files. 
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IRIS DMC MetaData Aggregator 
Station summary (1 time span) 


Legend: HQi^© 


Network 

YM :: TAIwan Integrated GEodynamic Research (TAIGERI:: YM Network Map :: DOI 

Station 

39 :: LDEO OBS 39 :: TAIGER :: 39 Station Map :: RESP :: SAC PZs :: XML 

Latitude 

21.599700 

Longitude 

121.364100 

Elevation 

-2508 

Start 

2008/05/11 (132)00:00:00 

End 

2009/05/30 (150) 23:59:59 

Epoch 

2008/05/11 (132) 21:25:51 - 2009/05/30 (150) 21:20:11 

Instrument 

LDEO OBS Sensor Mk2/LDEO_ncw dummy at 40/125 

Channels (Hz) 

Location BHl (40)0, BH2 (40)0, EHZ (40)0 

Instrument 

LDEO DPG/LDEO_ncw dummy at 40/125 

Channels (Hz) 

Location BDH (40]Q 

MetaData Load 

2015/03/30(089) 14:25:20 


Virtual network affiliations: 

Name Description Primary DC Secondary DC 

OBSIP U.S. National Ocean Bottom Seismograph Instruments OBSIP IRIS DMC 

■OBSIP LDEO OBSIP Experiments from LDEO OBSIP IRIS DMC 

OBSIP LDEO OBSIP Experiments from LDEO OBSIP IRIS DMC 

■UNRESTRICTED All unrestricted stations, generated via cron IRIS DMC IRIS DMC 

No data available in real-time systems 

Archive data availability - Make a batch request for data fbreg fast) - ( data access overviewi 

Earliest Latest 

□ 2(K)8/05/ll (132)22:25:51 2(X)9A)5/30(150)21:20:11 

Query archive for data day availability: 

2008 2009 


Figure 16. Detailecd Information on OBS Selectecd from 

Global Map. Source: [20]. 

C. REQUESTING DATA FROM IRIS 

There are multiple ways to request data from IRIS. The metho(d used in 
this thesis made use of the link in Figure 16, “Make a batch request for data 
(breq_fast).” The breq_fast request form, displayed in Figure 17, is self- 
explanatory with required items that must be entered starred. The “label” can be 
anything the user wants, which is the name of the file sent by the Data 
Management Center (DMC) to the entered email. The DMC is the part of IRIS 
that stores and makes available the seismic data [21]. After the “Start Query” 
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button is pressed, Figure 18 depicts the final webpage opened. For this thesis, 
“Full SEED” was the format used. After the form is submitted, IRIS sends an 
email with a link to download the requested SEED file directly. The link is 
available for 24 hours. “Full SEED” is used because of the various calibration 
files that Full SEED makes available with IRIS software. 


summaries 

by station 

by network 

by timesehes 

virtual nets 

breq.fast I 

channels 

stations 

responses 

temp networks 

assembled 

events comments 


BREQ_FAST Request Form 


virtual network 

network ym 


latitude and longitude 
NORTH 


station 39 
location 
channel 

data start time* 2008 Jun Q 13 Q ZE 000000 

data end time* 2008 Jun Q 13 Q ^ 2359 S 9 


sample rate >= and <= 
flags like 
comments like 
sensor type like 
site like 

data quality Best Q 


channel parameters 


WEST 

u 

EAST 





SOUTH 


Clear 

elevation 

depth 

>= 

>= 

<= 

<= 

azimuth 

dip 

>= 

>= 

<= 

<= 


name*: Jeremy Hankins 
institution*: Naval Postgraduate Schoo' 

Street address*: 1 university Cir., Monterey, 

email*: I 
phone*: 
label*: 

media: Electronic (FTP) Q 

alternate media: oat omm Q 

alternate media: olt Q 

DMC archived waveforms or 
query over ^ metadata/dataless/RESP** 

start Query 

•These fields arc required. 

♦*A query on metadata may result in a request file for unavailable data, but returns much faster. 


Figure 17. The BREQ_FAST Request Form Used to Get SEED Data 

from IRIS. Source: [22]. 
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Notice 


If this query takes more than 5 minutes, you may wish to hit the STOP button on this browser window, further restrict your query, and resubmit the 
query. 

Querying network like 'YM' for ycar=2008 
Found waveform data for 1 stations and 4 channels. 


Make Map of Stations 


SuDmit request below to receive 


.J full SEED 
datalcss SEED 
mini SEED 
sync file 
RESPfilc 

email request to myself 


.NAME Jeremy Hankins 

.INST Naval Postgraduate School 

.MAIL I University Cir., Monterey, CA 93943 

.EMAIL jhankinslnps.edu 

.PHONE 9045044708 

.PAX yOUR_FAX 

.MEDIA: FTP 

.ALTERNATE MEDIA: DAT 

.ALTERNATE MEDIA: DLT3 

.LABEL YM_39_13JUN08 

.FROM SO 

.QUALITY B 

.END 


39 

YM 

2008 

06 

13 

00 

00 

00.0000 

2008 

06 

13 

23 

59 

59.9999 

1 

BDH 

fView 

in. 

MPA) 

39 

YM 

2008 

06 

13 

00 

00 

00.0000 

2008 

06 

13 

23 

59 

59.9999 

1 

BHl 

(View 

Ail 

MPA) 

39 

YM 

2008 

06 

13 

00 

00 

00.0000 

2008 

06 

13 

23 

59 

59.9999 

1 

BH2 

(View 

in 

MDA1 

39 

YM 

2008 

06 

13 

00 

00 

00.0000 

2008 

06 

13 

23 

59 

59.9999 

1 

BHZ 

(View 

in 

MDA) 


Figure 18. Final Submission Page to Get SEED Data from IRIS. 

Source: [23]. 


D. DATA ANALYSIS 

It is at this point where a Mac computer (or another Unix-based machine) 
is required to process the data. As pointed out in Scobo’s thesis [1], Unix is 
required to run IRIS’s available software. A Mac computer proves useful because 
of the UNIX availability in the “Terminal” application. With an Internet search of 
Unix commands combined with the manuals that are available on IRIS’s website, 
a novice Unix user can at least use the installed programs from IRIS with 
practice. Installing and compiling Unix code, though, is not a trivial endeavor 
without previous Unix experience, and it is recommended to get assistance 
installing and compiling from someone with skills in Unix before proceeding with 
the data analysis portion. The IRIS website conveniently has instruction manuals 
for all of its software with examples that greatly assist novice operators. 
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Rdseed v5.3.1 


Rdseed reads the SEED files and outputs the raw data in SAC format, 
which is a binary format that can be read by the SAC software or MATLAB. In 
addition, RESP files and PZ files for all available channels can be output by 
rdseed for use in the calibration process. Rdseed proved effective and somewhat 
user-friendly after installed. Figure 19 displays a terminal window with typical 
rdseed commands. It is recommended to have an interactive manual open during 
the first few times using this software. The user can choose the dates and times 
to limit the SAC file’s size. A user could also download a multiple-day SEED file 
and then use rdseed to chop the SEED file into the desired times for SAC file 
output. The rdseed depicted in Figure 19 is in the prompt mode, which walks the 
user through the inputs to get a SAC file, a RESP file, and a PZ file (if desired). 


phs-imac:rdseedv5.3.1 phimac7$ ./rdseed 
« IRIS SEED Reader, Release 5.3.1 » 

Input File (/dev/nrste) or 'Quit’ to Exit: /users/phiniac7/desktop/YH_39/YH_39_DEC_08/YM_39_OE 
C_13_14.104940. seed 
Output File (stdout) 

Volume # [(1)-N) 

Options (acCsSpRtde] 

Summary file (None) 

Station List (ALL) 

Channel List (ALL) 

Network List (ALL) 

Loc Ids (ALL ("—" for spaces]) : 

Output Format [(1=SAC), 2=AH, 3=CSS, 4=mini seed, 5=seed, 6=SAC ASCII, 7=SEGY, 8=Siraple ASCIK 

SLIST), 9=Simple ASCII(TSPAIR)1 : 1 

Output file names include endtime? (Y/(N)]Y 

Output poles & zeroes ? [Y/(N))N 

Check Reversal [(0=No), l=Dip.Azimuth, 2=Gain, 3=Both]: 0 

Select Data Type ((E=Everything), O-Oata of Undetermined State, M=Merged data, R=Raw waveform 
Data, 0=QC'd data] :E 

Start Time(s) YYYY,ODD,SS.FFFF : 2008,348,00:00:00 
End Time(s) YYYY,ODD,SS.FFFF : 2008,348,23:59:59.9999 
Sample Buffer Length (20000000): 

Extract Responses [Y/(N)] : Y 

Output data format will be sac.binary. 

INFO: sac variable EVDP will be in KILOMETERS 

Writing YM.39..BDH, 3456000 samples (binary), starting 2008,348 00:00:00.0080 UT 

Writing YH.39..BH1, 3456000 samples (binary), starting 2008,348 00:00:00.0080 UT 

Writing YM.39..BH2, 3456000 samples (binary), starting 2008,348 00:00:00.0080 UT 

Warning... Azimuth/Dip Reversal found 39.BH2, Data inversion was not selected 

Writing YH.39..BHZ, 3456000 samples (binary), starting 2008,348 00:00:00.0080 UT 

Writing RESPONSE file: RESP.YH.39..BOH 
Writing RESPONSE file: RESP.YH.39..8H1 

.BH2 
.8HZ 


Writing RESPONSE file: RESP.YH.39.. 
Writing RESPONSE file: RESP.YH.39. 


Input File (/dev/nrst0) or 'Quit' to Exit: [] 


Figure 19. Example Rdseed Usage in Prompt Mode. 
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2 . 


EVALRESP 


For this thesis, EVALRESP was used to produce frequency amplitude 
phase (FAP) files. The FAP file is an ASCII file with three columns (frequency, 
amplitude, and phase) which can be used to apply the sensor response and 
obtain calibrated data. EVALRESP is built into the SAC program’s “transfer” 
command, but if the user wants to obtain a FAP file, the EVALRESP software 
must be run separately. According to the SAC Manual [24], the advantage of 
using a FAP with the “transfer” function “is that one can include additional stages 
of the instrument response and/or control more explicitly the frequency range 
over which the correction is applied.” A frequency range is chosen when building 
a FAP file within EVALRESP. 

Figure 20 depicts a sample EVALRESP command to generate a FAP file 
from the YM 39 OBS used in this thesis. Because channels 1, 2, and Z were 
used in data analysis, a FAP file for all three channels had to be generated in the 
method shown in Figure 20 for channel Z. The EVALRESP manual [25] provides 
detailed information on the commands in Figure 20. In this figure, “39” is the 
station, “BFIZ” is the channel, 2008 is the year, 348 is the Julian date, “0 20’ lists 
the frequency range sought after, “2049” is the requested number of samples in 
the given frequency range, and “-f” precedes the file location. The RESP file 
called, in this case “RESP.YM.39..BFIZ,” is the response file generated by the 
rdseed program, “-s” is the type of spacing requested (“lin” for linear spacing is 
chosen here), and “-r” is the response type which is set to “FAP” [25]. 

phs-imac:src phimac7S ./evalresp 39 BHZ 2008 348 0 20 2049 -f /Users/phiniac7/desktop/YH_39/YH_39 
_DEC_08/RESP.YH.39..BHZ -s tin -r 'fap' 
phs-imac:src phittiac7$ | 

Figure 20. Sample EVALRESP Input for FAP File. 

In addition to a FAP file, EVALRESP can also output a complex spectra 
(CS) file and an amplitude phase (AP) file. The CS file is a three column file 
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containing the frequency, real part of the amplitude, and imaginary part of the 
amplitude, whereas the AP file is similar to the FAP file except that frequency is 
not included. 

3. SAC 

SAC software, with the built-in transfer function, produced the calibrated 
data for this thesis. Figure 21 demonstrates the commands used for this thesis. 
Again, it is recommended to open an interactive software manual, available on 
the IRIS website, when using SAC. After reading-in the SAC data file with the 
read (r) command, any linear trends and offsets were removed with a remove 
trend command (rtr). The digital Fourier transform (DFT) that is used with a 
transfer function can pad the beginning and end of a time series with zeros so 
that SAC can achieve a power of 2 number of points. The “taper” command was 
applied to the data sequence after removing the trend to reduce sidelobes in the 
spectral data. The default setting for the taper command is a Planning taper with 
a width of .05 times the entire length of data used [24]. Therefore, the amplitude 
of 10% of the data sequence is reduced by the taper and was discarded from 
analysis. 

The SAC program has multiple built-in seismic analysis tools. 
Unfortunately, the SAC manual does not go into much detail discussing the 
theory behind how these analysis tools operate. The “transfer” command built 
into the SAC program “performs deconvolution to remove an instrument 
response and convolution to apply another instrument response.” [24] The 
transfer command can be achieved using either a RESP file, a FAP file, or a PZ 
file. During a transfer with FAP file option, if a frequency is encountered below 
the lowest frequency in the FAP file, the algorithm uses the lowest frequency 
given for the correction. Similarly, if the frequency encountered is above the max 
FAP file frequency, the algorithm uses the highest frequency given. The rdseed 
program produces the PZ and RESP files. The EVALRESP program, as 
discussed earlier, produces the FAP file. The four numbers following “freq” in 
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Figure 21 are the frequency limits, f1 through f4. According to the SAC Manual 
[24], the frequency limits should be such that f1 < f2 < f3 < f4. A high-pass taper 
is applied between f1 and f2, no taper (unity) is applied between f2 and f3, and a 
low-pass taper is applied between f3 and f4. To avoid ringing, [24] recommends 
that f2=2*f1 and f4>=2*f3. When using SAC, the file that is read-in will remain 
there unless the program is commanded to do otherwise. The INICM command 
deletes any data in memory and was used prior to any new file to be read in [24]. 


SEISMIC ANALYSIS CODE [11/11/2013 (Version 101.6a)] 

Copyright 1995 Regents of the University of California 

SAC> r /users/phiiiiac7/desktop/7D_J65A/NOV2011/NOV15/2011.319.00.00.00.0195_2011.319.2 
3.59.59.0195.7D.J65A..BHZ.M.SAC 
SAC> rtr 
SAC> taper 

SAC> trans from fap s /users/phifflac7/desktop/FAP.70.J65A..BHZ to vel freq .01 .02 10 
20 

Extrapolating below lowest FAPFILE frequency 
Extrapolating above highest FAPFILE frequency 
Station (J65A ), Channel (BHZ ) 

SAC> write NOV15_2011_7D_J65A_BHZ_fap01_to_25_calibrated_01_02_10_20.sac 
SAC> INICM 

Figure 21. SAC Example of Transfer Command Using FAP File to 

Calibrate Data. 


4. RSAC.m 

The rsac.m program is a publicly available MATLAB program that reads a 
SAC file and outputs a three column MATLAB file that is easily useable for 
analysis. The first rsac.m column in the output is time, the second column is 
amplitude, and the third column contains the SAC header information and was 
not used in this thesis. 

5. MATLAB Data Analysis 

MATLAB was used for data analysis on the calibrated SAC files after the 

data was read using rsac.m. Three functions (see Appendix A) were written that 

take as inputs a maximum desired amplitude, the sampling frequency of the 

OBS, and the data from the three orthogonal seismometer channels from an 

OBS. The first function entitled “vibration_analysis.m” finds the peaks of the 
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signal above the maximum amplitude and checks the times between these 
peaks. It then looks for the longest period of time between peaks. If a quiescent 
period (a period of time during which the signal remains below the maximum 
amplitude) is found in the vertical channel that is at least an hour long, a power 
spectral density (PSD) using the pwelch function in MATLAB is conducted on the 
quiescent period. The PSD displays all three orthogonal channels. The data 
within the taper at the beginning and end of the input time series is not included 
in the quiescent period and is removed prior to a PSD being conducted. If there 
is not at least an hour of quiescent time, a different MATLAB function was used. 
If the signal was too high and no hour-long quiescent periods exist, the function 
high_amp_vibration_analysis.m was used. This program finds the three-hour 
time period with the lowest mean amplitude and outputs the same information as 
the vibration_analysis.m program. If the 24-hour sample was low amplitude and 
did not cross the maximum amplitude line, the low_amp_vibration_analysis.m 
was conducted, which outputs a PSD for the entire sample. 

There are also three different types of histograms produced by these 
programs. The first is a histogram of the amplitudes throughout a 24-hour period. 
Both non-logarithmic and logarithmic y-axes were required to adequately display 
all of the occurrences. The next two histograms involved the desired maximum 
amplitude threshold which is entered into the functions. The first is a histogram of 
width of pulses that exceed the maximum amplitude threshold, and the second is 
a histogram of lengths of time between these pulses that exceed the maximum 
amplitude threshold. Both of these histograms are completed logarithmically and 
non-logarithmically to ensure that outliers are displayed. 

The maximum amplitude threshold histograms were completed because 
vibration amplitudes above a certain level would likely result in a noise floor too 
high to obtain satisfactory performance for the bottom mounted sensor. However, 
since the performance of a future bottom mounted sensor is unknown, the 
maximum amplitude threshold was chosen somewhat arbitrarily to be 5x10'^ m/s 
based on observations of high amplitude activity. This threshold proved to be a 
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good indicator of high amplitude, seismic activity (non-quiescent), and the 
threshold was thus used to compare the velocities measured at all of the OBS 
locations. Only YM 39 had standard deviations for 24-hour periods that regularly 
exceeded this threshold, and the other three OBSs, with the exception of a 
couple of days, remained well below. 

E. CHALLENGES FACED WITH CALIBRATION 

The largest challenge faced for this thesis was obtaining calibrated data 
before analysis could be conducted in MATLAB. [8] describes the LDEO OBS as 
able to effectively measure from .01 to 20 Hz. The LDEO OBS’s response curve, 
illustrated in Figure 22, changes by 2 orders of magnitude from .01 to 1 Hz. Figure 
23 exhibits the SIO OBS’s response curve that has nearly a flat response from .01 
to 20 Hz. Because of the LDEO response curve, getting the calibrated data 
becomes more complicated than just dividing the raw counts by the sensitivity. 


YM.39..BHZ 

2008-05-llT21:25:51 to 2009-05-30T21:20:11 



Figure 22. LDEO OBS Response Curve for YM 39 Channel Z. 

Source: [26]. 
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7D.J65A..BHZ 

2011-10-17722:31:00 to 2012-07-16707:21:00 



Figure 23. SIO OBS Response Curve from Channel Z of 7D J65A. 

Source: [27]. 

Due to lack of familiarity with SAC software, and because the SAC manual 
does not go into detail about the process behind calibration, calibration was 
initially conducted with MATLAB using the CS file. This proved to be more 
complicated than necessary, especially after conversation with a seismologist 
who stated that the SAC software was quite able to complete the calibration 
correction sufficiently. There were four ways attempted to complete the 
calibration that involved the IRIS SAC software and the time series Webserver 
(TSWS). 

Figure 24 displays the TSWS webpage. The TSWS allows a user to pick 
the network, station, channel, date(s), time, and multiple calibration/filtration 
options. Depending on the selection, it will immediately produce a plot or 
downloadable data in a user-selected format. For rsac.m to work with the TSWS 
data, “SAC binary little-endian” must be selected. With the correct settings, the 
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TSWS provided a convenient way to double check time plots. Without correct 
settings, the TSWS produced some odd results that could be misleading. 

To illustrate, Figure 25 is a TSWS plot of a full day from YM 39, with 
frequency limits of .01-.02-10-20 Hz, the instrument correction applied, and 
“remove mean” selected. Figure 26 is a MATLAB plot produced with FAR 
calibrated data of the same time period and settings. Both plots appear similar. 


Web Services IRISWS Timeseries Docs v. 1 Builder 


URL Builder, timeseries vl 


Service interface URL Builder Help Revisions 

Use this form to buM a URL to the timeseries wreb service. Notice that as you edit the form, the link is automaticaly updated. 


O Usage 


Network; 

70 

Remove mean; 

o 

Station; 

J6SA 

Low-Pass Filter: 

■ 

Location; 


High-Pass Fitter: 


Channel: 

BHZ 

Band-Pass Filter; 


Start Time; 

201MM5T12;00;00 ■ 

Differentiate: 


End Time: 

201MM5T12;01;21 ■ 

Integrate: 


Correction: 

Perform Instrument Correction ; 

Envelope: 


Frequency 

O .0V.02-1B-20 

Taper; 


Limits: 






Decimate (samples per 


Units: 

Scale: 


sec): 

Output: 

SAC binary little-endian ; 


CHch the link: 

http://service-iris.edu/iriswsAimeseries/1/query?net-7D&sta-J65A4cha-BHZ&start-2011-11-15T12:00:00Aend-201M1-15T12:01:21&freqllmits-.0K02-10-204 

demean-true&output-sacbl&k)C“-&corTect«tfue 


This form primarily serves to illustrate how to create ws-timeseries URLs and includes a fixed sequence of signal processing options. The ws-timeseries 
servKe performs all sijpial processing options in the order they appear within the URL. this flenbility is thus not offered by the form; advanced users 
may wish to manually modify or create their own request URLs with the processing options in their preferred order. 


Figure 24. TSWS Webpage. Source: [28]. 
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YM.39.7?.BHZ. 3456000 samples, 40.0 SpS, 2008-06-13X00:00:00.002 



Figure 25. TSWS Plot of YM 39 on June 13, 2008, 24 Hours, 
Frequency Limits of .01-.02-10-20 Hz. Source: [28]. 


YM3913JUN 08 All Day 



Figure 26. MATLAB Plot of YM 39 on June 13, 2008, 24 Hours, 
Frequency Limits of .01-.02-10-20 Hz, FAP Calibrated. 
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However, if the frequency range is changed slightly to 0-.02-10-20 Hz, and 
all other settings remain unchanged, the TSWS produces the plot in Figure 27. 
Figure 28 displays the FAR calibrated file with the same settings. Although the 
plot has changed with apparent ringing, it still is recognizable as a seismogram 
with close to the same values as the previous FAR plot. 


YM.39.7P.BH2, 3456000 samples, 40,0 sps, 2008-06-13T00:00;00.002 



Figure 27. TSWS Riot of YM 39 on June 13, 2008, 24 Hours, 
Frequency Limits of 0-.02-10-20 Hz. Source: [28]. 
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YM 3913 JUN 08 All Day 



Figure 28. MATLAB Plot of YM 39 on June 13, 2008, 24 Hours, 

Frequency Limits of 0-.02-10-20 Hz, FAP Calibrated. 

Smaller samples of the same day were inspected with similar results. 
Figure 29 shows a TSWS output of 4080 samples, from 12:00:00 to 12:01:42 on 
13 June 2008. The settings were the same as for Figure 25, with frequency limits 
of .01-.02-10-20 Hz, instrument correction applied, and “remove mean” selected. 
Figure 30 displays what happens to this smaller data sample when “remove 
mean” is not selected with all other settings unchanged. 
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(XlO-S) 


YM.39.??.BHZ, 4080 samples, 40.0 SpS, 2008-06-13T12;00:00.001 



Figure 29. TSWS Plot of YM 39 on June 13, 2008, from 12:00:00 to 
12:01:42, Frequency Limits of .01-.02-10-20 FIz. Source: [28]. 


YM.39.??.BHZ. 4080 samples, 40 0 sps, 2008-06-13112 00:00.001 



Figure 30. TSWS Plot of YM 39 on June 13, 2008, from 12:00:00 to 
12:01:42, Frequency Limits of .01-.02-10-20 Hz, Remove Mean 
Deselected. Source: [28]. 


39 






Figure 31 is a FAR calibrated plot of the same time period without “remove 
trend” applied. This behavior of the TSWS led to the conclusion that other 
calibration methods were needed besides using the user-friendly TSWS. If the 
frequency limits are applied correctly, the instrument correction is applied, and 
“remove mean” is selected, the amplitude of TSWS plots can still differ from the 
FAR, RZ, and RESR calibrated plots, but the shape of the signal is nearly 
identical to the other calibration methods. 


YM 39 13 JUN 08 12:00:00-12:01;42 



Figure 31. MATLAB Riot of YM 39 on June 13, 2008, 12:00:00 to 
12:01:42, Frequency Limits of .01-.02-10-20 FIz, FAR 
Calibrated, Remove Trend Not Completed. 


The remaining SAC calibration methods produced very similar results if 
the frequency limits were set the same. It was initially thought that FAR produced 
superior results to both a RESR and RZ calibrated file, but it was later discovered 
that FAR, RESR, and RZ calibrated files are nearly identical. Figures 32 and 33 
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show four different methods to gain calibrated data for both OBSs. Figure 32 is 
from YM 39 on June 13, 2008, and Figure 33 is from 7D J65A on November 15, 
2011. Of note, all the plots are nearly identical, but the TSWS produces different 
amplitudes in Figure 32. 


YM 3913 JUN 08 1200:00-12:01:42 
RESP Calibrated 



YM 3913 JUN 08 1200:00-1201:42 



YM 3913 JUN 08 1200:00-1201:42 



YM 3913 JUN 08 1200:00-1201:42 



Figure 32. MATLAB Plot of Different Calibration Methods for 

YM 39 on June 13, 2008. 
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7D J65A 15 NOV 2011 12:00:00-12X)1:21 



TifT>e(s) .10* 


7D J65A15 NOV2011 12:00:00-1201:21 



7D J65A15NOV2011 12:00:00-1201:21 7D J65A15 NOV 2011 12:00:00-1201:21 




Figure 33. MATLAB Plot of Different Calibration Methods for 
7D J65A on November 15, 2011. 

Because of the similarity of the data from the three SAC calibration 
techniques and because of the work already completed with FAP files, the FAP 
calibration technique was used for the data in this thesis. Any of the methods 
available with the SAC transfer, either FAP, RESP, or PZ files, would have been 
sufficient for use. 

Some idiosyncrasies with the software proved to hold additional 
challenges. Each time a transfer command in SAC was used, the command line 
involved a “from” type of file (i.e. RESP, PZ, FAP) and an optional “to” desired 
units (e.g., velocity and displacement). Figure 21 is an example of transfer with a 
“to vel” (to velocity) command. The default waveform is displacement, so if a “to 
none” line is entered, the SAC manual [24] states that the output will default to 
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displacement. If no “to” line is entered (left blank), SAC will assume it is “to none” 
[24]. A transfer with a FAP file, though, does not necessarily follow the SAC 
manual, which led to an assumption that calibrated files were in units of velocity 
when the files output were actually acceleration. Figure 34 summarizes the SAC 
peculiarities with the transfer command that led to some analysis challenges. 



POLEZERO VEL DISP DISP 

FAP ACC VEL VEL 


Figure 34. Summary of SAC Transfer Output. 

F. SUMMARIZED DATA ANALYSIS STEPS 

The steps described in this chapter are summarized in Figure 35. 



Figure 35. Methodology and Data Analysis Summary. 
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V. RESULTS 


A. POWER SPECTRAL DENSITY 

1. YM39andYM38 

YM 39 was the most erratic of the OBSs sampled. It had the highest 
variation in standard deviation and the largest variation in percentage of time 
spent above max amplitude. Most samples did not have at least one hour of 
quiescent time and had to be analyzed with the high amplitude vibration analysis 
code. In most instances for YM 39, the PSDs taken from the quiescent periods 
displayed a decreasing amplitude with increasing frequency and no apparent 
narrowband frequencies as illustrated in Figure 36; however, the decibel levels 
differed from sample to sample. Discrete frequencies do show up in the PSDs of 
the entire 24-hour period, which is displayed in Figure 37. 


Power Spectral Density Estimate 
Three Hours with Lowest Mean Amplitude 



Figure 36. YM 39, PSD from August 1,2008, Quiescent Period (Three 

Flours with Lowest Mean Amplitude). 
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Power Spectral Density Estimate 
Entire Day of Data 



Figure 37. YM 39, August 1,2008, PSD for Entire Day 

If the PSD for channel Z is graphed on a logarithmic x and y axis, it appears 
almost as a straight line as in Figure 38. Most spectral densities for YM 39 
channel Z follow this pattern with only slight variations from the overall slope. The 
velocity time series plot responsible for this PSD, Figure 39, represents a typical 
appearance of the majority of the days for YM 39. 
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Power Spectral Density Estimate 



Figure 38. YM 39, August 1,2008, PSD with Logarithmic Axes. 



Figure 39. YM 39 Velocity Time Series, August 1,2008. 

For the low frequency range, shown in Figure 40, most of YM 39’s PSDs 
exhibited a steady decline with the same overall pattern from month to month. 
Flowever, even though the overall shape and slope remained consistent, the 
decibel level varied greatly from sample to sample and was far more inconsistent 
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than the other OBSs. The drop off in decibel levels close to 0 Hz is due to the 
bandpass filter, which starts with a taper at .01 Hz. 


Power Spectral Density Estimate 
Three Hours with Lowest Mean Amplitude, 0 to 2 Hz 



Figure 40. YM 39, PSD, August 1,2008, 0 to 2 Hz. 

YM 38, however, had the lowest amplitudes of all four OBSs, especially in 
the vertical channel. Channel 1 and channel 2 exhibited much higher amplitudes 
than channel Z for all days sampled. This is shown by the fact that the decibel 
level of the vertical channel hovered around 10 to 20 dB less than the horizontal 
channels. The channel Z plot showed hints of discrete frequencies, especially 
around 4 Hz, 6 Hz, and 8 Hz, but the appearance of the discrete frequencies 
differed slightly from month to month like in Figure 41 where the discrete 
frequencies are located at approximately 4 Hz and 7 Hz. A visual comparison of 
plots between YM 38 and YM 39 on the same dates revealed no common 
discrete frequencies and different appearances of the time series plots. In 
particular, YM 38 was quieter than YM 39 by around 2 orders of magnitude. 
Figure 41 displays a typical YM 38 PSD for a quiescent period. Figure 42 is from 
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the same sample between 0 to 2 Hz. Most 24-hour samples from YM 38 never 
had any amplitude spikes above the max amplitude threshold. 


Power Spectral Density Estimate 
Longest Continuous Quiescent Period 



Figure 41. YM 38, PSD for Quiescent Period, August 1,2008. 


Power Spectral Density Estimate 
Longest Continuous Quiescent Period, 0 to 2 Hz 



Frequency (Hz) 


Figure 42. YM 38, PSD for Quiescent Period, August 1,2008, 0 to 2 

Hz. 
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The sample from August 1,2008 is an example of one of the few days when 
the amplitude crossed the maximum amplitude threshold. As displayed in Figure 
43, there are two large amplitude occurrences. Figure 44 is a close-up of the 
largest amplitude spike, and the amplitude increase does not appear to cause 
any ringing on the plot. 



Figure 43. YM 38, Velocity Time Series, August 1,2008. 


CIOM-iip Of High Amplitude, 1 AUG 2008, YM 38, Vertical Velocity (Channel Z) 



Figure 44. YM 38, August 1,2008, Close-Up of Flighest Amplitude 

Event in Velocity Time Series 
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Unlike in OBS YM 39, most of the high amplitudes in the horizontal 
channels do not show up in the vertical channel. Figure 45 is a depiction of 
channel 1 from YM 38 on August 1,2008. The higher amplitudes in the beginning 
of the day do not translate into a high response in channel Z, but the two high- 
amplitude occurrences after the 3x10"^ second mark and 4x10"^ second mark do 
appear in the vertical channel. 



. EntItvOayof Data with 0«sir*d Maximum Amplitud* Displayed, Horizontal V«lodty(ChanrMl1) 

• - 


Tun# lit 


Figure 45. YM 38 Velocity Time Series, Channel 1, August 1,2008. 

Although YM 38 and YM 39 were located in the Luzon Strait in close 
proximity, the two OBSs displayed vastly different results. YM 39 was not 
consistent, and YM 38 displayed much lower vertical vibrational disturbances. 

2. 7D J65A and 7D J73A 

Both OBSs west of the Strait of Juan de Fuca, 7D J65A and 7D J73A 
display similar and consistent results. 7D J65A exhibits the same consistent 
pattern almost every day where channel Z peaked at 8 to 9 FIz with varying 
decibel levels. Figures 46 and 47 depict the two 7D J65A plots with maximum 
and minimum decibel levels at around 9 Hz. Channel Z in Figure 46 peaks at 
around -130 dB and in Figure 47 peaks at around -150 dB. 
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Power Spectral Density Estimate 
Three Hours with Lowest Mean Amplitude 



Figure 46. 7D J65A, PSD of Quiescent Period (Three Hours with 
Lowest Mean Amplitude), February 15, 2012. 


Power Spectral Density Estimate 
Three Hours with Lowest Mean Amplitude 



Figure 47. 7D J65A, PSD of Quiescent Period (Three Hours with 
Lowest Mean Amplitude), March 1,2012. 
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Figure 48 is from 7D J73A, March 1, 2012, and shows an increase in 
decibels after 10 Hz. Even though a low-pass taper is applied at 10 Hz, 7D 
J73A’s decibel levels continue to increase after this taper is applied. Decibel 
levels for channel Z peak after the taper from as low as -160 dB to as high 
as -130 dB. This pattern of decibel level increase after the taper is present in 
nearly every sample. 


Power Spectral Density Estimate 
Three Hours with Lowest Mean Amplitude 



Figure 48. 7D J73A, PSD of Quiescent Period (Three Hours with 
Lowest Mean Amplitude), March 1,2012. 


From 0 to 2 Hz, 7D J65A and 7D J73A are similar, with both OBSs peaking 
at around 0.2 Hz and then displaying a steady, almost linear, decline after the 
peak. Figures 49 and 50 depict 0 to 2 Hz from both OBSs on the same day, 
March 1,2012. 
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Power Spectral Density Estimate 
Three Hours with Lowest Mean Amplitude, 0 to 2 Hz 



Figure 49. 7D J65A, PSD of Quiescent Period (Three Hours with 
Lowest Mean Amplitude), March 1,2012, 0 to 2 Hz. 


Power Spectral Density Estimate 
Three Hours with Lowest Mean Amplitude, 0 to 2 Hz 



Figure 50. 7D J73A, PSD of Quiescent Period (Three Hours with 
Lowest Mean Amplitude), March 1,2012, 0 to 2 Hz. 
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PSD Conclusion 


Figures 51 and 52 summarize the OBS PSD variability, focusing in on 1 
Hz and 5 Hz. The decibel levels for both 1 Hz and 5 Hz are graphed per sample 
for all samples taken. The results show fairly steady decibel levels for all OBSs at 
the frequencies except YM 39, which is the least consistent and noisiest of all the 
OBSs sampled. The decibel levels at 5 Hz are similar for YM 38 and the two 
OBSs located west of the Strait of Juan de Fuca. At 1 Hz, though, 7D J65A and 
7D J73A are consistently stronger. Tables listing all the values in Figures 51 and 
52 are located in Annex B. 


PSD Values at 1 Hz and S Hz, YM 39 and YM 38 

14-«ufv4C l.M-0} ]S-Jul4« i-AugOt IS-Auf-Oi t-Mp-Oi l-Oct-O* tS-Oct-W l-NpvM tS-N»v«a iS-Ow-W l-Mn-09 t-M«y 0* lVM*r.09 



Figure 51. Ouiescent PSD Values at 1 Hz and 5 Hz for 

YM 39 and YM 38. 


PSD Values at 1 Hz and 5 Hz, 70 J65A and 7D J73A 

l-Movll IS-Novll I-Om- 11 lS^-11 I-Ma- 12 IS-Jart-lZ 1-F«b-12 lS-rtb-12 1-M«r-12 lS-M«r-12 l-Apf-12 1$-Apr-12 1-Mtvl2 15 -M«y- 12 l-Jun-12 15-Jun*12 l-Kil-12 lS-Jur-12 
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-•-70 )65A dS At 1 Hz -#-70 J6SA dB At 5 Hz -•-70 i73A dB At 1 Hz 


-•-70J73AdBAtSHz 


Figure 52. Ouiescent PSD Values at 1 Hz and 5 Hz for 

7D J65A and J73A. 
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Although acceleration is not discussed in Results, acceleration can be 
obtained from the velocity data used in all of the figures by using Equation 21, 
where a is acceleration, co is angular frequency, and v is velocity. Appendix B 
contains rms acceleration for all vertical channels. 

a - icov ( 21 ) 

B. AMPLITUDE HISTOGRAMS 
1. 7D J65A and 7D J73A 

Amplitude histograms were generally viewed with a logarithmic y-axis so 
that more details were available. However, any velocity time series plot that was 
consistent and without multiple large amplitude spikes appeared Gaussian if 
viewed without a logarithmic axis. Most of the 7D J65A and 7D J73A amplitude 
histograms were Gaussian in appearance. Every amplitude histogram from these 
two OBSs viewed on a logarithmic scale looked similar to either Figure 53 or 
Figure 56. Figure 53 is from OBS 7D J73A, November 1, 2011. Figure 54 
displays the same histogram, but with a non-logarithmic y-axis. The channel Z 
time series plot represented by these histograms is shown in Figure 55. 
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Histogram of ChannsI Z Amplitudos 



Figure 53. 7D J73A, Amplitude Histogram, November 1,2011. 



-1.5 -1 -0.5 0 0£ 1 15 2 

VelocftyAmpIftudes (nrVs) .iq-* 


Figure 54. 7D J73A, Amplitude Histogram, November 1,2011, Non- 

Logarithmic. 
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Entire Day of Data with Desired Maximum Amplitude Displayed, Vertical Velocity (Channel Z) 
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Figure 55. 7D J73A, Velocity Time Series, November 1,2011. 

Figure 56 displays a narrower histogram in a logarithmic plot. Such narrow 
histograms were also common in many of the days from 7D J65A and 7D J73A. 
The histogram in Figure 56 naturally has a smaller number of amplitude events 
that surpass the amplitude threshold. Figures 57 and 58 display the same 
information in non-logarithmic plots, but Figure 58 is a close-up displaying the 
single, larger amplitude events that occurred throughout this sample day. Figure 
59 is the velocity time series plot of channel Z for July 1, 2012, represented by 
the histograms in Figures 56 through 58. Both OBSs are consistent enough that 
the wider and narrow histograms in Figure 53 and Figure 56 are representative of 
the range of standard deviation and shapes of all the data investigated. 
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Figure 56. 7D J73A, Amplitude Histogram, July 1,2012. 
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Figure 57. 7D J73A, Amplitude Histogram, July 1,2012, Non- 

Logarithmic. 
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Figure 58. 7D J73A, Amplitude Histogram, July 1,2012, 

Close-Up Non-Logarithmic. 
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Entire Day of Data with Desired Maximum Amplitude Displayed, Vertical Velocity (Channel Z) 
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Figure 59. 7D J73A, Velocity Time Series, July 1,2012, 

Channel Z. 


2. YM 39 and YM 38 

For the YM 39 and YM 38 OBSs, the non-logarithmic amplitude 
histograms miss many details otherwise available with a logarithmic plot. Many 
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time series plots for both YM 38 and YM 39 involve long, relatively quiet time 
periods followed by shorter periods of higher amplitude vibrations. Figure 60 from 
June 13, 2008, displays this type of time series plot. The amplitude histogram in 
Figure 61 is logarithmic, and the amplitude histogram in Figure 62 is not and 
misses useful specifics provided in Figure 61. Many of the histograms for YM 38 
and YM 39 have similar appearances, although the samples from YM 38 were up 
to two orders of magnitude smaller in standard deviation. 



Figure 60. YM39, Velocity Time Series, June 13, 2008, Channel Z. 
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Figure 61. YM 39, Amplitude Histogram, June 13, 2008. 



Figure 62. YM 39 Amplitude Histogram, Non-Logarithmic, 

June 13, 2008. 
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Multiple amplitude histograms from both YM 38 and YM 39 appear in a 
crown-type shape shown in Figure 63. Figure 64 is the time series plot for the 
same day, February 1, 2009. This histogram illustrates the common shape for 
the plots when there were calm periods interrupted by multiple, almost impulsive, 
amplitude spikes. Figures 65 and 66, from August 1, 2008, display the same 
behavior with YM 38 but at a much smaller amplitude. Initially this behavior was 
thought to be electronic clipping, but there are histogram examples with the same 
appearance but increased amplitude. The vibrations arrive and are symmetric 
with an equal positive and negative amplitude response. 


. Histogram of Qtannel Z Amplitudes 
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Figure 63. YM 39, Amplitude Flistogram, February 1,2009. 
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Figure 64. YM 39, Velocity Time Series, February 1,2009. 
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Figure 65. YM 38, Amplitude Flistogram, August 1,2008. 
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Figure 66. YM 38, Velocity Time Series, August 1,2008, Channel Z. 

3. Amplitude Histogram Conclusion 

The histograms for both the OBSs located west of the Strait of Juan de 
Fuca proved that the data these two OBSs collected was consistent with few 
outliers. This led to amplitude histograms with the appearance of a normal 
distribution. Many samples taken from the OBSs in the Luzon Strait displayed 
long, relatively quiet periods, which were abruptly interrupted by higher amplitude 
seismic readings. These particular events were responsible for the crown-type 
shapes in Figures 63 and 65. 

C. HISTOGRAMS OF HIGH AMPLITUDE PULSE WIDTHS AND 

HISTOGRAMS OF TIME BETWEEN HIGH AMPLITUDE PULSES 

Two types of histograms were based on a maximum entered amplitude 
threshold: a histogram of the width of pulses that exceeded the maximum 
amplitude and a histogram of times between these pulses. As mentioned in 
Chapter IV, the maximum amplitude threshold chosen for this thesis was 5x1 O’® 
m/s. Although this threshold was higher than the standard deviation of any OBS 
except for YM 39, all of the OBSs have occurrences where the amplitude 
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surpasses this line. This threshold proved to be a good indicator of high 
amplitude, non-quiescent activity. 

1. 7D J65A and 7D J73A 

Both OBSs west of the Strait of Juan de Fuca spent around 6% of the 
sampled time above the maximum amplitude threshold. For these OBSs, the 
maximum signal length above the maximum amplitude threshold was around 6 
seconds, but the majority of these high amplitude signals were 1 second or less 
in length. Figure 67 is a representative histogram of signals that exceeded the 
maximum amplitude threshold for both 7D J65A and 7D J73A. Figure 68 is the 
time series associated with this histogram. 


Histogram of Length of Continuous Time Signal Exceeding Maximum 
Channel Z 
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Figure 67. 7D J65A, November 15, 2011, Width of Signals 
Exceeding Maximum Amplitude Threshold 
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Figure 68. 7D J65A, November 15, 2011, 

Velocity Time Series 

The vast majority of times between signals that exceeded the maximum 
amplitude was 1 second or less. Figures 69 and 70 display time between events, 
and Figure 71 is a close up with a bin width of 1 second. The higher amplitude 
signals occur in groups, which explains the large number of occurrences of the 1 
second or less bin. 


Histogram of Length of Time Between Signals that Exceed Max Amplitude 
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Figure 69. 7D J65A, November 15, 2011, Time Between 
Signals That Exceed Maximum Amplitude 
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Histogram of Length of Time Between Signals that Exceed Max Amplitude 
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Figure 70. 7D J65A, November 15, 2011, Time Between Signals That 
Exceed Maximum Amplitude, Non-Logarithmic 
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Figure 71. 7D J65A, November 15, 2011, Time Between 
Signals That Exceed Maximum Amplitude 
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YM 39 AND YM 38 


OBS YM 39 spent an average of 52% of time above the maximum 
amplitude threshold over all the days sampled, although the percentage of time 
spent above the threshold varied significantly from 1% to 91% per day sampled. 
The sample from March 15, 2009, spent 54% of the time above the threshold, 
which is close to the overall YM 39 average of 52%. Figure 72 depicts the time 
series of channel Z from March 15, 2009. Figure 73 is the histogram of signal 
widths above the maximum amplitude threshold, and Figure 74 displays the 
times between these signals. Figure 75 is non-logarithmic, and displays the 
single occurrences where the largest time period between signals is almost five 
minutes. 



Figure 72. YM 39, March 15, 2009, Velocity Time Series 
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Figure 73. YM 39, March 15, 2009, Width of Signals Exceeding 

Maximum Amplitude Threshold 


Histogram of Length of Time Between Signals that Exceed Max Amplitude 

Channel Z 



Figure 74. YM 39, March 15, 2009, Time Between Signals That Exceed 

Maximum Amplitude Threshold 
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Histogram of Length of Time Between Signais that Exceed Max Ampiitude 
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Figure 75. YM 39, March 15, 2009, Time Between Signals 
That Exceed Maximum Amplitude Threshold 

3. Conclusion 

It was anticipated that the length of time between high amplitude signals 
would be short for the days with a high percentage of time spent above the 
maximum amplitude threshold. What was unexpected, however, was the fact that 
the majority of times between high amplitude signals were 1 second or less for 
an OBS that spent only 7% of its time above the maximum amplitude (Figure 68). 
This confirms that the high amplitude signals tend to occur in groups. The 
maximum time spent between signals on this same day was 16 minutes, but this 
only occurred once in a 24 hour sample. Table 2 summarizes the extremes of the 
OBSs. 
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Table 2. Maximum and Minimum Standard Deviations and Times 
Spent Above the Maximum Amplitude for All OBSs 


OBS 

Max standard deviation 
observed (i-im/s) 

Min standard deviation 
observed (nm/s) 

Max % of time spent 
above maxamplitude 

Min % of time spent 
above maxamplitude 

YM 39 

662 

16.1 

91 

1 

YM 38 

7.3 

0.23 

0.16 

0 

7DJ65A 

117 

7.2 

24 

0 

7DJ73A 

39.9 

7.6 

21 

0 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. METHODOLOGY 

The global view on the IRIS website proved useful in the selection of 
OBSs that were co-located in significant geographic chokepoints. IRIS also made 
it convenient to select the OBS directly, view the detailed information, and 
request the data directly from the OBS information webpage. The publicly 
available seismic processing software on the IRIS website was essential in 
obtaining calibrated data from the raw data. Manual calibration was attempted, 
but the IRIS software, which was designed by the seismology community, was 
quite capable once familiarity with the software was established. There were 
multiple paths to calibrated data, especially using the SAC transfer function, but 
the novice user will need to ensure that the SAC transfer output is in the units 
that he/she is expecting. A Mac computer, or a Unix-based operating system, is 
required to operate the IRIS software. The MATLAB rsac.m function was a 
convenient method to convert SAC data to a MATLAB format. 

B. RESULTS 

Out of the four OBSs, three proved consistent from month to month in 
their individual discrete frequencies and decibel levels, but there were no 
common discrete frequencies between OBSs that were located in the same 
geographic area. YM 39 in the Luzon Strait was erratic with large differences in 
decibel levels from sample to sample. The amplitude histograms for 7D J65A and 
7D J73A displayed what appeared to be normal distributions, with the exception 
of a couple of outliers with some seismic activity noted. YM 39 and YM 38 
displayed oddly shaped amplitude histograms that were the result of long 
quiescent periods interrupted by large, impulsive seismic activity. The histograms 
based on the maximum amplitude threshold showed that the higher amplitude 
vibrations occurred in groups, with the vast majority of counts occurring with a 
second or less between events. Even the histograms from samples where the 
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percentage of time spent above the maximum amplitude threshold was low 
showed results where the vast majority of time between signals was one second 
or less. 

The extreme differences between YM 38 and YM 39 lead to questions 
about environmental effects on YM 39 that are not affecting YM 38. Even though 
YM 38 exhibited some high-amplitude activity on occasion, it was far more quiet 
and consistent than YM 39. There was not much in the way of commonality 
between these two OBSs despite their close proximity in the same geographic 
area. 

7D J65A and 7D J73A had close values in decibel levels, standard 
deviation, and percentage of time spent above the maximum amplitude threshold 
in almost every sample taken on the same day, with the exception of a couple of 
outliers higher amplitude seismic events. It would appear that west of the Strait of 
Juan de Fuca is a fairly homogeneous location with regard to environmental 
factors on the ocean floor. 

C. RECOMMENDATIONS FOR FUTURE RESEARCH 

1. Transfer Functions to Reduce Long Period Noise 

Pressure from infragravity waves on the surface of the ocean can 
penetrate to the ocean floor causing deformation of the ocean floor [5]. According 
to [5], these pressure variations are the source of long period seismic noise, from 
around 30 seconds to 100 or more seconds (1 to 3.3 mHz). These pressure 
variations are also picked up by the DPG. Therefore, a transfer function can be 
developed that can get rid of the seismic noise caused by the surface waves. 
This could also possibly be used to take noise out of any future bottom mounted 
sensor. 

2. Environmental Effects in the Luzon Strait 

The extreme differences between the two OBSs in the Luzon Strait 
present questions about environmental effects on YM 39 that are not affecting 
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YM 38. The Luzon Strait is known for its large internal wave activity, and 
investigating a correlation between this activity with OBS results would be an 
interesting study. A detailed study of currents and tides in the local region 
combined with OBS results could also prove useful. 

3. Correlation of Data with Other OBS in Luzon Strait, Possible 
Mechanical Failure of YM 39 

Other OBSs in the neighborhood of YM 38 and YM 39 could be inspected 
to correlate the results obtained in this thesis. IRIS data problem reports were 
searched, but no OBSs in the YM network were listed. Currently, then, there are 
no indications of a mechanical failure of YM 39, but additional data could be 
explored to both confirm the results and rule out any mechanical failure of YM 39. 

4. Electronic Noise Floor of OBSs 

The electronic noise floors of the OBSs were not found during the 
research for this thesis, and it is possible that the amplitude histograms with a 
Gaussian appearance were the result of OBS electronic noise. Thus, the OBSs 
could be limited by the electronic noise floor during quiescent periods. No 
instrument noise models were found that provide the rough order of magnitude of 
the electronic noise floor, and any future study should explore these questions 
further. In particular, it would be valuable to have confidence in the ability of OBS 
sensors to accurately measure ocean accelerations down to 10 ® m/s^ without 
running into the electronic noise floor. 

5. OBS Requests 

The OBSIP website [4] contains information for researchers to request 
OBSIP instrument usage, and the OBSIP allows their OBS equipment to be 
available to other government research or educational institutions [4]. Therefore, 
LLNL or NPS could request usage of OBSs in other areas of interest. 
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APPENDIX A. MATLAB CODE 


A. VIBRATION ANALYSIS CODE 

function vibration_analysis(rsac_channel_Z,max_amp,Fs,rsac_channeM ,rsac_channel_2) 
% OBS, ocean bottom seismometer, maximum ampiitude, BHZ, sampiing frequency 
% 

% Purpose: Takes three caiibrated channeis of OBS data, OBS sampiing 
% frequency, and maximum desired ampiitude, and finds the maximum 
% quiescent period (at ieast 1 hour in iength where the ampiitude 
% does not exceed maximum ampiitude). Function outputs an ampiitude (m/s) 

% vs time piot (for aii three channeis), the mean and standard deviation 
% of entire sampie ampiitude, power spectrai densities of iongest 
% quiescent period, a percentage of time spent above maximum desired 
% ampiitude, a histogram of ampiitudes, a histogram of iengths of signais 
% greater than max desired ampiitude, and a histogram of time periods 
% between signais that exceed maximum ampiitude. 

% 

% Output Variables: 

% std_z = standard deviation of channel z amplitude 
% percentage_time_above_max_amp = percentage of time where amplitude 
% is greater than max amplitude 
% 

% Input Variables: 

% rsac_channel_Z = vertical .sac file after it has been read through 
% rsac.m 

% max_amp = desired maximum amplitude (m/s) 

% Fs = sampling frequency of OBS (Hz) 

% rsac_channel_1 = channel 1 .sac file after it has been read through 
% rsac.m 

% rsac_channel_2 = channel 2 .sac file after it has been read through 
% rsac.m 
% 

% Local Variables: 

% cal_data = calibrated vertical channel amplitude (m/s) 

% time = time (seconds) 

% PKS = value of peaks greater than max amplitudel (m/s) 

% LOGS = vector of location of PKS (integer location) 

% loopjndex = binary file (index) for location of PKS 

% loopjndexjogical = logical version of loopjndex 
% cal_data1 = same as cal_data, but modified to zero-out peak 
% locations (m/s) 

% PKS2 = value of peaks in second peaks-finding function (m/s) 

% LOCS2 = vector of location of PKS2 
% loopJndexJogical2=logical version of modified loopjndex 
% cal_data2 = same as cal_data1, but modified to zero-out additional 
% peaks (m/s) 

% cal_data3 = same as cal_data2, but modified to contain longest 
% quiescent occurrence only (m/s) 

% largest_quiescentjndex = index locating the largest quiescent 
% period 

% channel1_caLdata = calibrated amplitude of channel 1 (m/s) 
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% channel2_cal_data = calibrated amplitude of channel 2 (m/s) 

% time_quiescent = time location of longest quiescent period (s) 

% N = FFT length 
% overlap = pwelch overlap value 

% Pxx = channel Z distribution of power per unit of freq, but in this 
% application units of [m/s/sqrt(Hz)] 

% F, Ft, F2 = Frequency (FIz) 

% Pxx1 = channel 1 distribution of power per unit of freq, but in this 
% application units of [m/s/sqrt(Hz)] 

% Pxx2 = channel 2 distribution of power per unit of freq, but in this 
% application units of [m/s/sqrt(Hz)] 

% index = index locating where cal_data is greater than the maximum 
% amplitude 

% above_max_amp_time = time where cal_data is greater than maximum 
% amplitude (s) 

% outlier = logical index location where Z-amplitude is greater than 
% maximum amplitude 

% cal_data_update = same as cal_data, but only contains amplitudes 
% greater than maximum amplitude, (m/s) 

% W = vector with lengths of signals that exceed max amplitude 
% S = vector that contains different lengths of time between signals 
% that exceed max amplitude 

% 

% Functions Called: None. 

% 

% Filename: vibration_analysis.m 
% Written by: Jeremy R. Hankins 

cal_data=rsac_channel_Z(:,2); 

time=rsac_channel_Z(:,1); 

cal_data_qplot=cal_data; 

rsac1_qplot=rsac_channel_1(:,2); 

rsac2_qplot=rsac_channel_2(:,2); 

% plot amplitude (m/s) vs. time (s) for all channels 

figure 

plot(time,cal_data) 
hold on 

plot([time(1) time(length(time))],[max_amp max_amp],'r-') 
hold on 

plot([time(1) time(length(time))],[-max_amp -max_amp],'r-') 
ylabel('Velocity (m/s)','Fontsize',12) 
xlabel('Time (s)','Fontsize',12) 

title('Entire Day of Data with Desired Maximum Amplitude Displayed, Vertical Velocity (Channel 
Z)','Fontsize',16) 

legend('Calibrated amplitude','+/- max amplitude') 
figure 

plot(time,rsac_channeM (:,2)) 
hold on 

plot([time(1) time(length(time))],[max_amp max_amp],'r-') 
hold on 

plot([time(1) time(length(time))],[-max_amp -max_amp],'r-') 
ylabel('Velocity (m/s)','Fontsize',12) 
xlabel('Time (s)','Fontsize',12) 
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title('Entire Day of Data with Desired Maximum Ampiitude Displayed, Horizontal Velocity 

(Channel 1)yFontsize',16) 

legend('Calibrated amplitude','+/- max amplitude') 

figure 

plot(time,rsac_channel_2(:,2)) 
hold on 

plot([time(1) time(length(time))],[max_amp max_amp],'r-') 
hold on 

plot([time(1) time(length(time))],[-max_amp -max_amp],'r-') 
ylabel('Velocity (m/s)','Fontsize',12) 
xlabel('Time (s)','Fontsize',12) 

title('Entire Day of Data with Desired Maximum Amplitude Displayed, Horizontal Velocity 

(Channel 2)','Fontsize',16) 

legend('Calibrated amplitude','+/- max amplitude') 

rsac_channel_1=rsac_channel_1(:,2); 

rsac_channel_2=rsac_channel_2(:,2); 

% find the peaks above the maximum amplitude 

[PKS,LOCS]=findpeaks(abs(cal_data),'MinPeakHeight',max_amp); 

% for-loop to create index locating where time between peaks is at least an 
% hour 

for i=1 :length(LOCS); 
if i==1 

if time(LOCS(i))-time(1 )>=(60*60); 

loop_index(1 :LOCS(i))=0; 
else loop_index(1 :LOCS(i))=1; 
end 

elseif i>1 && i<length(LOCS); 

if time(LOCS(i))-time(LOCS(i-1 ))>=(60*60) 
loop_index(LOCS(i-1):LOCS(i))=0; 
else loop_index(LOCS(i-1):LOCS(i))=1; 
end 

elseif i==length(LOCS); 

if length(time)-time(LOCS(i))>=(60*60) 
loop_index(LOCS(i):length(time))=0; 
else loop_index(LOCS(i):length(time))=1; 
end 
end 
end 

loopJndex_logical=logical(loopJndex); 
cal_data1 =cal_data; 
cal_data1(loopjndex_logical)=0; 

% get rid of remaining peaks (which are located directly next to previous 
% peak locations) 
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[PKS2,LOCS2]=findpeaks(abs(cal_data1),'MinPeakHeight',max_amp); 
for i=1 :length(LOCS2) 

loop_index(LOCS2(i)-60*Fs:LOCS2(i)+60*Fs)=1; %taking away 1 minute on either side of 
remaining peaks 
end 

ioop_index_logical2=logical(loop_index); 
cai_data2=cai_data1; 
cai_data2(loop_index_logical2)=0; 

% find largest quiescent period 
cal_data3=cal_data2; 

largest_quiescent_index=ones(length(cal_data2),1); 

X=diff(LOCS); 

X=[LOCS(1)-1; X; length(cal_data2)-LOCS(length(LOCS))]; % need to add beginning length 
l=find(X==max(X)); 


if l==1 && length(time(round(length(time)*.05):LOCS(l)))>=60*60*Fs 
largest_quiescent_index(round(length(time)*.05):LOCS(l))=0; 
elseif l==1 && length(time(round(length(time)*.05):LOCS(l)))<60*60*Fs 
l=find(X==max(X(2:length(X)))); 
end 

if l==length(X) && length(time(LOCS(l-1):length(cal_data2)-round(.05*length(time))))>=60*60*Fs 
largest_quiescent_index(LOCS(l-1):length(cal_data2)-round(.05*length(time)))=0; 
elseif l==length(X) && length(time(LOCS(l-1):length(cal_data2)- 
round(.05*length(time))))<60*60*Fs 
l=find(X==max(X(2:length(X)-1))); 
end 

if l>1 && l<length(X); 

largest_quiescent_index(LOCS(l-1):LOCS(l))=0; 

end 

% Ensure X is at least an hour 
Y = X(I); 
if Y<=60*60*Fs 

error='quiscent period not an hour use high amp vibration analysis' 
elseif Y>60*60*Fs 
display('hour of quiescent found') 
end 

largest_quiescent_index=logical(largest_quiescentJndex); 

cal_data3(largest_quiescent_index)=0; 

channel1_cal_data=rsac_channel_1; %channel 1 
channel1_cal_data(largest_quiescentjndex)=0; 

channel2_cal_data=rsac_channel_2; %channel2 
channel2_cal_data(largest_quiescent_index)=0; 
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time_quiescent=time(~largest_quiescentjndex); 

% plot longest quiescent period 
figure 

plot(tlme_qulescent,cal_data_qplot(~largest_qulescent_lndex)) 
hold on 

plot([tlme_qulescent(1) tlme_qulescent(length(tlme_qulescent))],[max_amp max_amp],'r-') 
hold on 

plot([tlme_qulescent(1) tlme_qulescent(length(tlme_qulescent))],[-max_amp -max_amp],'r-') 
hold on 

plot(tlme_qulescent,rsac1_qplot(~largest_qulescent_lndex),'g') 
hold on 

plot(tlme_qulescent,rsac2_qplot(~largest_qulescent_lndex),'m') 
ylabel('Veloclty (m/s)','Fontslze',12) 
xlabel(Tlme (s)','Fontslze',12) 

tltle('Longest Continuous Quiescent Period (Based on Channel Z)','Fontslze',16) 
legend('Callbrated amplitude, channel Z','+/- max amplitude','+/- max amplitude','Calibrated 
amplitude, channel 1','Calibrated amplitude, channel 2') 

% Power Spectral Density of quiescent period 


N=2^12; 

overlap=N/2; 

figure 

[Pxx,F]=pwelch(cal_data3(~largest_qulescentJndex),hammlng(N),overlap,N,Fs); 

[Pxx1 ,F1]=pwelch(channel1_cal_data(~largest_qulescentJndex),hammlng(N),overlap,N,Fs); 
[Pxx2,F2]=pwelch(channel2_cal_data(~largest_qulescent_lndex),hammlng(N),overlap,N,Fs); 
plot(F,10*log10(Pxx),'r') 
hold on 

plot(F1,10*log10(Pxx1),'b') 
hold on 

plot(F2,10*log10(Pxx2),'m') 

ylabel('dB re m/s/$$\sqrt{Hz}$$','Interpreter','Latex','Fontsize', 14) 
xlabel('Frequency (Hz)','Fontslze',12) 

tltle({'Power Spectral Density Estimate','Longest Continuous Quiescent Period'},'Fontsize', 16) 
grid on 

legend('Channel Z','Channel 1','Channel 2') 
figure 

plot(F,10*log10(Pxx),'r') 
hold on 

plot(F1,10*log10(Pxx1 ),'b') 
hold on 

plot(F2,10*log10(Pxx2),'m') 

ylabel('dB re m/s/$$\sqrt{Flz}$$','Interpreter','Latex','Fontsize', 14) 
xlabel('Frequency (FIz)','Fontsize', 12) 
xllm([0 2]) 

tltle({'Power Spectral Density Estimate','Longest Continuous Quiescent Period, 0 to 2 
Hz'}, 'Fontsize', 16) 
grid on 

legend('Channel Z','Channel 1','Channel 2') 


length_z=length(cal_data); 
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taperz=floor(.05*length_z); 


length_1 =length(rsac_channeM); 
taperl =floor(.05*length_1); 

length_2=length(rsac_channel_2); 

taper2=floor(.05*length_2); 

cal_data=cal_data(taperz:length_z-taperz); 
rsac_channel_1 =rsac_channel_1 (taperl :length_1 -taperl); 
rsac_channel_2=rsac_channel_2(taper2:length_2-taper2); 

time=time(taperz:length_z-taperz); 

% Find percentage of time above maximum ampiitude 

index=abs(cal_data)>=max_amp; 

above_max_amp_time=time(index); 

percentage_time_above_max_amp=(length(above_max_amp_time)/length(time))*100 

% Create square wave where amplitude surpasses maximum amplitude for use in 
% histogram. 

cal_data_update=cal_data; 
outlier=abs(cal_data_update) > max_amp; 

cal_data_update(~outlier)=0; % Sets entire signal below max amp to zero 
cal_data_update(outlier)=sign(cal_data_update(outlier)); % This sets all 
% magnitudes above max amp to 1 (builds a square wave) 

W=pulsewidth(abs(cal_data_update),Fs); % stores times of pulses above max amp 

% W = pulsewidth(X) returns a vector, W, containing the time differences 
% between the mid-reference level instants of the initial and final transitions 
% of each positive-polarity pulse in the bilevel waveform, X. 

[S,INITCROSS]=pulsesep(abs(cal_data_update),Fs); 

S=[INITCROSS(1); S]; 

% S = pulsesep(X) returns the differences, S, between the mid-reference level 
%instants of the final negative-going transitions of every positive-polarity 
%pulse and the next positive-going transition. X is a bilevel waveform. 

% [S,INITCROSS] = pulsesep(...) returns the mid-reference level instants, 

% INITCROSS, of the first positive-polarity transitions. 

figure 

histogram(W) 

set(gca,'Yscale','log') % this gives you average time of pulses above max amp 
title({'Flistogram of Length of Continuous Time Signal Exceeding Maximum','Channel 
Z'},'Fontsize',16) 
xlabel('Time (s)','Fontsize',12) 

figure 

histogram(W) 
ylim([0 50]) 
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xlabel('Time (s)','Fontsize',12) 

title({'Histogram of Length of Continuous Time Signai Exceeding Maximum','Channei 
Z'},'Fontsize',16) 

figure 

histogram(S) 
set(gca,'Yscaie','iog') 
xiabei('Time (s)','Fontsize',12) 

titie({'Flistogram of Length of Time Between Signais that Exceed Max Amplitude','Channei 
Z'},'Fontsize',16) 

figure 

histogram(S) 
ylim([0 50]) 

xlabel('Time (s)','Fontsize',12) 

titleff'Histogram of Length of Time Between Signals that Exceed Max Amplitude','Channel 
Z'},'Fontsize',16) 

figure 

histogram(cal_data,400) 
set(gca,'Yscale','log') 
xlabel('Velocity Amplitude (m/s)') 
ylabel('Number of Occurrences') 
title('Flistogram of Channel Z Amplitudes') 

figure 

histogram(cal_data,400) 
xlabel('Velocity Amplitude (m/s)') 
ylabel('Number of Occurrences maxed at 20') 
ylim([0 20]) 

title('Flistogram of Channel Z Amplitudes') 
std_z=std(cal_data) 

B. HIGH AMPLITUDE VIBRATION ANALYSIS CODE 

function 

high_amp_vibration_analysis(rsac_channel_Z,max_amp,Fs,rsac_channeM,rsac_channel_2) 
% OBS, ocean bottom seismometer, maximum amplitude, BFIZ, sampling frequency 
% 

% Purpose: Run this function if the sample day too noisy for 
% vibration_analysis.m (not at least 1 hour in length where the amplitude 
% does not exceed maximum amplitude). Takes three calibrated channels of 
% OBS data, OBS sampling frequency, and maximum desired amplitude, and 
% finds the three hours with lowest mean amplitude. Function outputs an 
% amplitude (m/s) vs time plot (for all three channels), the mean and 
% standard deviation of entire sample amplitude, power spectral densities 
% of three-hour with lowest mean, a percentage of time spent above maximum 
% desired amplitude, a histogram of amplitudes, a histogram of lengths of 
% signals greater than max desired amplitude, and a histogram of time 
% periods between signals that exceed maximum amplitude. 

% 

% Sample needs to be at least 15 hours long for this function to work 
% properly. 
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% 

% Output Variables: 

% std_z = standard deviation of channel z amplitude 
% percentage_time_above_max_amp = percentage of time where amplitude 
% is greater than max amplitude 
% 

% Input Variables: 

% rsac_channel_Z = vertical .sac file after it has been read through 
% rsac.m 

% max_amp = desired maximum amplitude (m/s) 

% Fs = sampling frequency of OBS (Hz) 

% rsac_channeM = channel 1 .sac file after it has been read through 
% rsac.m 

% rsac_channel_2 = channel 2 .sac file after it has been read through 
% rsac.m 
% 

% Local Variables: 

% cal_data = calibrated vertical channel amplitude (m/s) 

% time = time (seconds) 

% length_z = length of cal_data (vector length), channel Z 
% taperz = amount to subtract from beginning^nd of sample to get rid 
% of taper, channel Z 

% length_1 = vector length of channel 1 amplitude vector 
% taperf = amount to subtract from beginning/end of sample to get rid 
% of taper, channel 1 

% length_2 = vector length of channel 2 amplitude vector 
% taper2 = amount to subtract from beginning/end of sample to get rid 
% of taper, channel 2 

% rsac_channeM = amplitude vector channel 1 

% rsac_channel_2 = amplitdue vector channel 2 

% second_3_hours_cal_data to twentyone_3_hours_cal_data = three-hour 
% blocks of time used to find three-hour time period with lowest mean 
% ampitidue 

% MEAN = vector containing mean values of three-hour blocks of time 

% Y = used in first if/elseif statement. If size of sample much less 

% than 24 hours, Y variable helps assign correct area to apply taper 
% X = used in if/elseif statement to choose correct three-hour block 
% and assign variables to PSD 

% N = FFT length 

% overlap = pwelch overlap value 

% Pxx = channel Z distribution of power per unit of freq, but in this 
% application units of [m/s/sqrt(Hz)] 

% F, F1, F2 = Frequency (Hz) 

% Pxx1 = channel 1 distribution of power per unit of freq, but in this 
% application units of [m/s/sqrt(Hz)] 

% Pxx2 = channel 2 distribution of power per unit of freq, but in this 
% application units of [m/s/sqrt(Hz)] 

% index = index locating where cal_data is greater than the maximum 

% amplitude 

% above_max_amp_time = time where cal_data is greater than maximum 
% amplitude (s) 

% outlier = logical index location where Z-amplitude is greater than 
% maximum amplitude 

% cal_data_update = same as cal_data, but only contains amplitudes 
% greater than maximum amplitude, (m/s) 
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% W = vector with lengths of signals that exceed max amplitude 
% S = vector that contains different lengths of time between signals 
% that exceed max amplitude 

% 

% Functions Called: None. 

% 

% Filename: high_amp_vibration_analysis.m 

% Written by: Jeremy R. Flankins 

cal_data=rsac_channel_Z(:,2); 

time=rsac_channel_Z(:,1); 

length_z=length(cal_data); 

taperz=floor(.05*length_z); 

length_1 =length(rsac_channel_1 (:,2)); 
taperl =floor(.05*length_1); 

length_2=length(rsac_channel_2(:,2)); 

taper2=floor(.05*length_2); 

rsac_channel_1=rsac_channel_1(:,2); 

rsac_channel_2=rsac_channel_2(:,2); 

figure 

plot(time,cal_data) 
hold on 

plot([time(1) time(length(time))],[max_amp max_amp],'r-') 
hold on 

plot([time(1) time(length(time))],[-max_amp -max_amp],'r-') 
ylabel('Velocity (m/s)','Fontsize',12) 
xlabel(Time (s)','Fontsize',12) 

title('Entire Day of Data with Desired Maximum Amplitude Displayed, Vertical Velocity (Channel 
Z)','Fontsize',16) 

legend('Calibrated amplitude','+/- max amplitude') 
figure 

plot(time,rsac_channel_1) 
hold on 

plot([time(1) time(length(time))],[max_amp max_amp],'r-') 
hold on 

plot([time(1) time(length(time))],[-max_amp -max_amp],'r-') 
ylabel('Velocity (m/s)','Fontsize',12) 
xlabel('Time (s)','Fontsize',12) 

title('Entire Day of Data with Desired Maximum Amplitude Displayed, Florizontal Velocity 

(Channel 1)','Fontsize',16) 

legend('Calibrated amplitude','+/- max amplitude') 

figure 

plot(time,rsac_channel_2) 
hold on 

plot([time(1) time(length(time))],[max_amp max_amp],'r-') 
hold on 

plot([time(1) time(length(time))],[-max_amp -max_amp],'r-') 
ylabel('Velocity (m/s)','Fontsize',12) 
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xlabel('Time (s)','Fontsize',12) 

title('Entire Day of Data with Desired Maximum Ampiitude Dispiayed, Horizontai Velocity 

(Channel 2)yFontsize',16) 

legend('Calibrated amplitude','+/- max amplitude') 

% Divide sample into 3 hour time blocks, exiuding taper at beginning and 
% end of sample. 3 hour block that includes taper could be less/more 
% than three hours. Amount less/more depends on length of sample. 

% if/elseif statement covers files as short as 15 hours 

second_3_hours_cal_data=cal_data(taperz:4*60*60*Fs); 

third_3_hours_cal_data=cal_data(2*60*60*Fs:5*60*60*Fs); 

fourth_3_hours_cal_data=cal_data(3*60*60*Fs:6*60*60*Fs); 

fifth_3_hours_cal_data=cal_data(4*60*60*Fs:7*60*60*Fs); 

sixth_3_hours_cal_data=cal_data(5*60*60*Fs:8*60*60*Fs); 

seventh_3_hours_cal_data=cal_data(6*60*60*Fs:9*60*60*Fs); 

eighth_3_hours_cal_data=cal_data(7*60*60*Fs:10*60*60*Fs); 

nine_3_hours_cal_data=cal_data(8*60*60*Fs:11 *60*60*Fs); 

ten_3_hours_cal_data=cal_data(9*60*60*Fs:12*60*60*Fs); 

eleven_3_hours_cal_data=cal_data(10*60*60*Fs:13*60*60*Fs); 

twelve_3_hours_cal_data=cal_data(11 *60*60*Fs:14*60*60*Fs); 

if length(cal_data)/(Fs*60*60)>=23; 

Y=20; 

thirteen_3_hours_cal_data=cal_data(12*60*60*Fs:15*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(13*60*60*Fs:16*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(15*60*60*Fs:18*60*60*Fs); 
seventeen_3_hours_cal_data=cal_data(16*60*60*Fs:19*60*60*Fs); 
eighteen_3_hours_cal_data=cal_data(17*60*60*Fs:20*60*60*Fs); 
nineteen_3_hours_cal_data=cal_data(18*60*60*Fs:21 *60*60*Fs); 
twenty_3_hours_cal_data=cal_data(19*60*60*Fs:22*60*60*Fs); 
twentyone_3_hours_cal_data=cal_data(20*60*60*Fs:length_z-taperz); 
elseif length(cal_data)/(Fs*60*60)>=22 && length(cal_data)/(Fs*60*60) < 23; 

Y=19; 

thirteen_3_hours_cal_data=cal_data(12*60*60*Fs:15*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(13*60*60*Fs:16*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(15*60*60*Fs:18*60*60*Fs); 
seventeen_3_hours_cal_data=cal_data(16*60*60*Fs:19*60*60*Fs); 
eighteen_3_hours_cal_data=cal_data(17*60*60*Fs:20*60*60*Fs); 
nineteen_3_hours_cal_data=cal_data(18*60*60*Fs:21 *60*60*Fs); 
twenty_3_hours_cal_data=cal_data(19*60*60*Fs:length_z-taperz); 
twentyone_3_hours_cal_data=1 El 0; 

elseif length(cal_data)/(Fs*60*60)>=21 && length(cal_data)/(Fs*60*60) < 22; 

Y=18; 

thirteen_3_hours_cal_data=cal_data(12*60*60*Fs:15*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(13*60*60*Fs:16*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(15*60*60*Fs:18*60*60*Fs); 
seventeen_3_hours_cal_data=cal_data(16*60*60*Fs:19*60*60*Fs); 
eighteen_3_hours_cal_data=cal_data(17*60*60*Fs:20*60*60*Fs); 
nineteen_3_hours_cal_data=cal_data(18*60*60*Fs:length_z-taperz); 
twenty_3_hours_cal_data=1 El 0; 
twentyone_3_hours_cal_data=1 El 0; 
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elseif length(cal_data)/(Fs*60*60)>=20 && length(cal_data)/(Fs*60*60) < 21; 
Y=17; 

thirteen_3_hours_cal_data=cal_data(12*60*60*Fs:15*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(13*60*60*Fs:16*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(15*60*60*Fs:18*60*60*Fs); 
seventeen_3_hours_cal_data=cal_data(16*60*60*Fs:19*60*60*Fs); 
eighteen_3_hours_cal_data=cal_data(17*60*60*Fs:length_z-taperz); 
nineteen_3_hours_cal_data=1 El 0; 
twenty_3_hours_cal_data=1 El 0; 
twentyone_3_hours_cal_data=1 El 0; 

elseif length(cal_data)/(Fs*60*60)>=19 && length(cal_data)/(Fs*60*60) < 20; 
Y=16; 

thirteen_3_hours_cal_data=cal_data(12*60*60*Fs:15*60*60*Fs); 

fourteen_3_hours_cal_data=cal_data(13*60*60*Fs:16*60*60*Fs); 

fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 

sixteen_3_hours_cal_data=cal_data(15*60*60*Fs:18*60*60*Fs); 

seventeen_3_hours_cal_data=cal_data(16*60*60*Fs:length_z-taperz); 

eighteen_3_hours_cal_data=1 El 0; 

nineteen_3_hours_cal_data=1 El 0; 

twenty_3_hours_cal_data=1 El 0; 

twentyone_3_hours_cal_data=1 El 0; 

elseif length(cal_data)/(Fs*60*60)>=18 && length(cal_data)/(Fs*60*60) < 19; 
Y=15; 

thirteen_3_hours_cal_data=cal_data(12*60*60*Fs:15*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(13*60*60*Fs:16*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(15*60*60*Fs:length_z-taperz); 
seventeen_3_hours_cal_data=1 El 0; 
eighteen_3_hours_cal_data=1 El 0; 
nineteen_3_hours_cal_data=1 El 0; 
twenty_3_hours_cal_data=1 El 0; 
twentyone_3_hours_cal_data=1 El 0; 

elseif length(cal_data)/(Fs*60*60)>=17 && length(cal_data)/(Fs*60*60) < 18; 
Y=14; 

thirteen_3_hours_cal_data=cal_data(12*60*60*Fs:15*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(13*60*60*Fs:16*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:length_z-taperz); 
sixteen_3_hours_cal_data=1 El 0; 
seventeen_3_hours_cal_data=1 El 0; 
eighteen_3_hours_cal_data=1 El 0; 
nineteen_3_hours_cal_data=1 El 0; 
twenty_3_hours_cal_data=1 El 0; 
twentyone_3_hours_cal_data=1 El 0; 

elseif length(cal_data)/(Fs*60*60)>=16 && length(cal_data)/(Fs*60*60) < 17; 
Y=13; 

thirteen_3_hours_cal_data=cal_data(12*60*60*Fs:15*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(13*60*60*Fs:length_z-taperz); 
fifteen_3_hours_cal_data=1 El 0; 
sixteen_3_hours_cal_data=1 El 0; 
seventeen_3_hours_cal_data=1 El 0; 
eighteen_3_hours_cal_data=1 El 0; 
nineteen_3_hours_cal_data=1 El 0; 
twenty_3_hours_cal_data=1 El 0; 
twentyone_3_hours_cal_data=1 El 0; 
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elseif length(cal_data)/(Fs*60*60)>=15 && length(cal_data)/(Fs*60*60) < 16; 

Y=12; 

thirteen_3_hours_cal_data=cal_data(12*60*60*Fs:length_z-taperz); 
fourteen_3_hours_cal_data=1 E10; 
fifteen_3_hours_cal_data=1 El 0; 
sixteen_3_hours_cal_data=1 El 0; 
seventeen_3_hours_cal_data=1 El 0; 
eighteen_3_hours_cal_data=1 El 0; 
nineteen_3_hours_cal_data=1 El 0; 
twenty_3_hours_cal_data=1 El 0; 
twentyone_3_hours_cal_data=1 El 0; 
end 

MEAN=[mean(abs(second_3_hours_cal_data)),... 

mean(abs(third_3_hours_cal_data)),mean(abs(fourth_3_hours_cal_data)),... 

mean(abs(fifth_3_hours_cal_data)),mean(abs(sixth_3_hours_cal_data)),... 

mean(abs(seventh_3_hours_cal_clata)),mean(abs(eighth_3_hours_cal_data)),... 

mean(abs(nine_3_hours_cal_data)),mean(abs(ten_3_hours_cal_data)),... 

mean(abs(eleven_3_hours_cal_data)),mean(abs(twelve_3_hours_cal_data)),... 

mean(abs(thirteen_3_hours_cal_data)),mean(abs(fourteen_3_hours_cal_data)),... 

mean(abs(fifteen_3_hours_cal_data)),mean(abs(sixteen_3_hours_cal_data)),... 

mean(abs(seventeen_3_hours_cal_data)),mean(abs(eighteen_3_hours_cal_data)),. 

mean(abs(nineteen_3_hours_cal_data)),mean(abs(twenty_3_hours_cal_data)),... 

mean(abs(twentyone_3_hours_cal_data))]; 

X=find(MEAN==min(MEAN)) 

% if/elseif statments to plot correct psd 


if X==1 

psd_z=second_3_hours_cal_data; 
psd_1 =rsac_channel_1 (taperl :4*60*60*Fs); 
psd_2=rsac_channel_2(taper2:4*60*60*Fs); 
elseif X==2 

psd_z=third_3_hours_cal_data; 
psd_1 =rsac_channel_1 (2*60*60*Fs:5*60*60*Fs); 
psd_2=rsac_channel_2(2*60*60*Fs:5*60*60*Fs); 
elseif X==3 

psd_z=fourth_3_hours_cal_data; 
psd_1 =rsac_channel_1 (3*60*60*Fs:6*60*60*Fs); 
psd_2=rsac_channel_2(3*60*60*Fs:6*60*60*Fs); 
elseif X==4 

psd_z=fifth_3_hours_cal_data; 
psd_1 =rsac_channel_1 (4*60*60*Fs:7*60*60*Fs); 
psd_2=rsac_channel_2(4*60*60*Fs:7*60*60*Fs); 
elseif X==5 

psd_z=sixth_3_hours_cal_data; 
psd_1 =rsac_channel_1 (5*60*60*Fs:8*60*60*Fs); 
psd_2=rsac_channel_2(5*60*60*Fs:8*60*60*Fs); 
elseif X==6 

psd_z=seventh_3_hours_cal_data; 
psd_1 =rsac_channeL1 (6*60*60*Fs:9*60*60*Fs); 
psd_2=rsac_channeL2(6*60*60*Fs:9*60*60*Fs); 
elseif X==7 

psd_z=eighth_3_hours_cal_data; 


88 


psd_1 =rsac_channel_1 (7*60*60*Fs:10*60*60*Fs); 
psd_2=rsac_channel_2(7*60*60*Fs:10*60*60*Fs); 
elseif X==8 

psd_z=n i n e_3_h o u rs_cal_d ata; 
psd_1 =rsac_channel_1 (8*60*60*Fs:11 *60*60*Fs); 
psd_2=rsac_channeL2(8*60*60*Fs:11*60*60*Fs); 
elseif X==9 

psd_z=ten_3_hours_cal_data; 
psd_1=rsac_channeL1 (9*60*60*Fs:12*60*60*Fs); 
psd_2=rsac_channeL2(9*60*60*Fs:12*60*60*Fs); 
elseif X==10 

psd_z=eleven_3_hours_cal_data; 
psd_1 =rsac_channeL1 (10*60*60*Fs:13*60*60*Fs); 
psd_2=rsac_channeL2(10*60*60*Fs:13*60*60*Fs); 
elseif X==11 

psd_z=twelve_3_hours_cal_data; 
psd_1=rsac_channel_1 (11*60*60*Fs:14*60*60*Fs); 
psd_2=rsac_channel_2(11*60*60*Fs:14*60*60*Fs); 
elseif X==12 

psd_z=thirteen_3_hours_cal_data; 
if Y==12 

psd_1 =rsac_channeL1 (12*60*60*Fs:length_1 -taperf); 
psd_2=rsac_channel_2(12*60*60*Fs;length_2-taper2); 
else 

psd_1 =rsac_channel_1 (12*60*60*Fs:15*60*60*Fs); 
psd_2=rsac_channel_2(12*60*60*Fs:15*60*60*Fs); 
end 

elseif X==13 

psd_z=fourteen_3_hours_cal_data; 
if Y==13 

psd_1 =rsac_channel_1 (13*60*60*Fs:length_1 -taperf); 
psd_2=rsac_channel_2(13*60*60*Fs:length_2-taper2); 
else 

psd_1 =rsac_channeL1 (13*60*60*Fs:16*60*60*Fs); 
psd_2=rsac_channeL2(13*60*60*Fs:16*60*60*Fs); 
end 

elseif X==14 

psd_z=fifteen_3_hours_cal_data; 
if Y==14 

psd_1 =rsac_channel_1 (14*60*60*Fs:length_1 -taperf); 
psd_2=rsac_channel_2(f 4*60*60*Fs:length_2-taper2); 
else 

psd_f =rsac_channel_f (f 4*60*60*Fs:f 7*60*60*Fs); 
psd_2=rsac_channel_2(f 4*60*60*Fs:f 7*60*60*Fs); 
end 

elseif X==f 5 

psd_z=sixteen_3_hours_cal_data; 
if Y==f5 

psd_f =rsac_channel_f (f 5*60*60*Fs:length_f -taperf); 
psd_2=rsac_channel_2(f 5*60*60*Fs;length_2-taper2); 
else 

psd_f =rsac_channeLf (f 5*60*60*Fs;f 8*60*60*Fs); 
psd_2=rsac_channeL2(f 5*60*60*Fs;f 8*60*60*Fs); 
end 

elseif X==f6 
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psd_z=seventeen_3_hours_cal_data; 
if Y==16 

psd_1 =rsac_channeL1 (16*60*60*Fs:length_1 -taperl); 
psd_2=rsac_channel_2(16*60*60*Fs;length_2-taper2); 
else 

psd_1 =rsac_channeL1 (16*60*60*Fs:19*60*60*Fs); 
psd_2=rsac_channeL2(16*60*60*Fs:19*60*60*Fs); 
end 

elseif X==17 

psd_z=eighteen_3_hours_cal_data; 
if Y==17 

psd_1 =rsac_channel_1 (17*60*60*Fs:length_1 -taperl); 
psd_2=rsac_channel_2(17*60*60*Fs:length_2-taper2); 
else 

psd_1 =rsac_channeL1 (17*60*60*Fs:20*60*60*Fs); 
psd_2=rsac_channeL2(17*60*60*Fs:20*60*60*Fs); 
end 

elseif X==18 

psd_z=nineteen_3_hours_cal_data; 
if Y==18 

psd_1 =rsac_channel_1 (18*60*60*Fs:length_1 -taperl); 
psd_2=rsac_channel_2(18*60*60*Fs;length_2-taper2); 
else 

psd_1 =rsac_channel_1 (18*60*60*Fs;21 *60*60*Fs); 
psd_2=rsac_channel_2(18*60*60*Fs:21 *60*60*Fs); 
end 

elseif X==19 

psd_z=twenty_3_hours_cal_data; 
if Y==19; 

psd_1 =rsac_channel_1 (19*60*60*Fs:length_1 -taperl); 
psd_2=rsac_channel_2(19*60*60*Fs:length_2-taper2); 
else 

psd_1 =rsac_channeL1 (19*60*60*Fs:22*60*60*Fs); 
psd_2=rsac_channeL2(19*60*60*Fs:22*60*60*Fs); 
end 

elseif X==20 

psd_z=twentyone_3_hours_cal_data; 

psd_1 =rsac_channeM (20*60*60*Fs:length_1 -taperl); 

psd_2=rsac_channel_2(20*60*60*Fs:length_2-taper2); 

end 


N=2^12; 

overlap=N/2; 

figure 

[Pxx,F]=pwelch(psd_z,hamming(N),overlap,N,Fs); 

[Pxx1 ,F1]=pwelch(psd_1 ,hamming(N),overlap,N,Fs); 
[Pxx2,F2]=pwelch(psd_2,hamming(N),overlap,N,Fs); 
plot(F,10*log10(Pxx),'r') 
hold on 

plot(F1,10*log10(Pxx1),'b') 
hold on 

plot(F2,10*log10(Pxx2),'m') 

ylabel('dB re m/s/$$\sqrt{Flz}$$','Interpreter','Latex','Fontsize', 14) 
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xlabel('Frequency (Hz)','Fontsize',12) 
xlim([0 2]) 

title({'Power Spectral Density Estimate','Three Flours with Lowest Mean Amplitude, 0 to 2 

Flz'},'Fontsize',16) 

grid on 

legend('Channel Z','Channel 1','Channel 2') 
figure 

plot(F,10*log10(Pxx),'r') 
hold on 

plot(F1,10*log10(Pxx1),'b') 
hold on 

plot(F2,10*log10(Pxx2),'m') 

ylabel('dB re m/s/$$\sqrt{Hz}$$','Interpreter','Latex','Fontsize', 14) 
xlabel('Frequency (Flz)','Fontsize',12) 

title({'Power Spectral Density Estimate','Three Hours with Lowest Mean Amplitude'}, 'Fontsize',1 6) 
grid on 

legend('Channel Z','Channel 1','Channel 2') 

% get rid of taper for plot and histograms 
cal_data=cal_data(taperz:length_z-taperz); 
rsac_channel_1 =rsac_channel_1 (taperl :length_1 -taperl); 
rsac_channel_2=rsac_channel_2(taper2:length_2-taper2); 
time=time(taperz:length_z-taperz); 

% psd for entire day of data 
figure 

[Pxx,F]=pwelch(cal_data,hamming(N),overlap, N,Fs); 

[Pxxl ,F1]=pwelch(rsac_channel_1,hamming(N),overlap,N,Fs); 
[Pxx2,F2]=pwelch(rsac_channel_2,hamming(N),overlap,N,Fs); 
plot(F,10*log10(Pxx),'r') 
hold on 

plot(F1,10*log10(Pxx1),'b') 
hold on 

plot(F2,10*log10(Pxx2),'m') 

ylabel('dB re m/s/$$\sqrt{Hz}$$','Interpreter','Latex','Fontsize', 14) 
xlabel('Frequency (Hz)','Fontsize',12) 

title({'Power Spectral Density Estimate','Entire Day of Data'}, 'Fontsize', 16) 
grid on 

legend('Channel Z','Channel 1','Channel 2') 

index=abs(cal_data)>=max_amp; 

above_max_amp_time=time(index); 

percentage_time_above_max_amp=(length(above_max_amp_time)/length(time))*100 

cal_data_update=cal_data; 
outlier=abs(cal_data_update) > max_amp; 

cal_data_update(~outlier)=0; % Sets entire signal below max amp to zero 

cal_data_update(outlier)=sign(cal_data_update(outlier)); % This sets all magnitudes above max 
amp to 1 (builds a square wave) 

W=pulsewidth(abs(cal_data_update),Fs); % stores times of pulses above max amp 

[S,INITCROSS]=pulsesep(abs(cal_data_update),Fs); 

S=[INITCROSS(1); Sj; 
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figure 

histogram(W) 

set(gca,'Yscale','log') % this gives you average time of puises above max amp 
titie({'Histogram of Length of Continuous Time Signal Exceeding Maximum','Channel 
Z'},'Fontsize',16) 
xlabel('Time (s)','Fontsize',12) 

figure 

histogram(W) 
ylim([0 50]) 

xlabel('Time (s)','Fontsize',12) 

title({'Histogram of Length of Continuous Time Signal Exceeding Maximum','Channel 
Z'},'Fontsize',16) 

figure 

histogram(S) 
set(gca,'Yscale','log') 
xlabel('Time (s)','Fontsize',12) 

title({'Histogram of Length of Time Between Signals that Exceed Max Amplitude','Channel 
Z'},'Fontsize',16) 

figure 

histogram(S) 
ylim([0 50]) 

xlabel('Time (s)','Fontsize',12) 

title({'Histogram of Length of Time Between Signals that Exceed Max Amplitude','Channel 
Z'},'Fontsize',16) 

figure 

histogram(cal_data,400) 
set(gca,'Yscale','log') 
xlabel('Velocity Amplitudes (m/s)') 
ylabel('Number of Occurrences') 
title('Histogram of Channel Z Amplitudes') 

figure 

histogram(cal_data,400) 
xlabel('Velocity Amplitudes (m/s)') 
ylabel('Number of Occurrences maxed at 20') 
ylim([0 20]) 

title('Histogram of Channel Z Amplitudes') 
std_z=std(cal_data) 
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c. 


LOW AMPLITUDE VIBRATION ANALYSIS CODE 


function 

low_amp_vibration_analysis(rsac_channel_Z,max_amp,Fs,rsac_channel_1,rsac_channel_2) 
% OBS, ocean bottom seismometer, maximum ampiitude, BHZ, sampiing frequency 
% 

% Purpose: Run this function if the sampie day does not exceed max 
% ampiitude with vibration_anaiysis.m (does not cross max ampiitude at 
% aii). Takes three caiibrated channeis of OBS data, OBS sampiing 
% frequency, and maximum desired ampiitude. Function outputs an 
% ampiitude (m/s) vs time piot (for aii three channeis), the mean and 
% standard deviation of entire sampie ampiitude, power spectrai densities 
% of entire day, and a histogram of ampiitudes. 

% 

% Sampie needs to be at ieast 15 hours iong for this function to work 
% properiy. 

% 

% Output Variables: 

% avg_z = mean of absolute value of channel z amplitude 
% std_z = standard deviation of absolute value of channel z amplitude 
% percentage_time_above_max_amp = percentage of time above max amp 
% (should be zero) 

% 

% Input Variables: 

% rsac_channel_Z = vertical .sac file after it has been read through 
% rsac.m 

% max_amp = desired maximum amplitude (m/s) 

% Fs = sampling frequency of OBS (Hz) 

% rsac_channel_1 = channel 1 .sac file after it has been read through 
% rsac.m 

% rsac_channel_2 = channel 2 .sac file after it has been read through 
% rsac.m 
% 

% Local Variables: 

% cal_data = calibrated vertical channel amplitude (m/s) 

% time = time (seconds) 

% length_z = length of cal_data (vector length), channel Z 
% taperz = amount to subtract from beginning/end of sample to get rid 
% of taper, channel Z 

% length_1 = vector length of channel 1 amplitude vector 
% tapert = amount to subtract from beginning/end of sample to get rid 
% of taper, channel 1 

% length_2 = vector length of channel 2 amplitude vector 
% taper2 = amount to subtract from beginning/end of sample to get rid 
% of taper, channel 2 

% rsac_channeM = amplitude vector channel 1 
% rsac_channel_2 = amplitdue vector channel 2 
% N = FFT length 
% overlap = pwelch overlap value 

% Pxx = channel Z distribution of power per unit of freq, but in this 
% application units of [m/s/sqrt(Hz)] 

% F, Ft, F2 = Frequency (Hz) 

% Pxx1 = channel 1 distribution of power per unit of freq, but in this 
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% application units of [m/s/sqrt(Hz)] 

% Pxx2 = channel 2 distribution of power per unit of freq, but in this 
% application units of [m/s/sqrt(Hz)] 

% index = index locating where cal_data is greater than the maximum 
% amplitude 

% above_max_amp_time = time where cal_data is greater than maximum 
% amplitude (s) 

% 

% Functions Called: None. 

% 

% Filename: low_amp_vibration_analysis.m 

% Written by: Jeremy R. Hankins 

cal_data=rsac_channel_Z(:,2); 

time=rsac_channel_Z(:,1); 

length_z=length(cal_data); 

taperz=floor(.05*length_z); 

length_1 =length(rsac_channeM (:,2)); 
taperl =floor(.05*length_1); 

length_2=length(rsac_channel_2(:,2)); 

taper2=floor(.05*length_2); 

rsac_channel_1=rsac_channel_1(:,2); 

rsac_channel_2=rsac_channel_2(:,2); 

figure 

plot(time,cal_data) 
hold on 

plot([time(1) time(length(time))],[max_amp max_amp],'r-') 
hold on 

plot([time(1) time(length(time))],[-max_amp -max_amp],'r-') 
ylabel('Velocity (m/s)','Fontsize',12) 
xlabel(Time (s)','Fontsize',12) 

title('Entire Day of Data with Desired Maximum Amplitude Displayed, Vertical Velocity (Channel 
Z)','Fontsize',16) 

legend('Calibrated amplitude','+/- max amplitude') 
figure 

plot(time,rsac_channeM) 
hold on 

plot([time(1) time(length(time))],[max_amp max_amp],'r-') 
hold on 

plot([time(1) time(length(time))],[-max_amp -max_amp],'r-') 
ylabel('Velocity (m/s)','Fontsize',12) 
xlabel('Time (s)','Fontsize',12) 

title('Entire Day of Data with Desired Maximum Amplitude Displayed, Horizontal Velocity 

(Channel 1)','Fontsize',16) 

legend('Calibrated amplitude','+/- max amplitude') 

figure 

plot(time,rsac_channel_2) 
hold on 
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plot([time(1) time(length(time))],[max_amp max_amp],'r-') 
hold on 

plot([time(1) time(length(time))],[-max_amp -max_amp],'r-') 
ylabel('Velocity (m/s)','Fontsize',12) 
xlabel('Time (s)','Fontsize',12) 

title('Entire Day of Data with Desired Maximum Amplitude Displayed, Horizontal Velocity 

(Channel 2)','Fontsize',16) 

legend('Calibrated amplitude','+/- max amplitude') 

% get rid of tapers 

cal_data=cal_data(taperz:length_z-taperz); 
rsac_channel_1 =rsac_channel_1 (tapert :length_1 -tapert); 
rsac_channel_2=rsac_channel_2(taper2:length_2-taper2); 
time=time(taperz:length_z-taperz); 


N=2^12; 

overlap=N/2; 

figure 

[Pxx,F]=pwelch(cal_data,hamming(N),overlap, N,Fs); 

[Pxx1 ,F1]=pwelch(rsac_channel_1,hamming(N),overlap,N,Fs); 
[Pxx2,F2]=pwelch(rsac_channel_2,hamming(N),overlap,N,Fs); 
plot(F,10*log10(Pxx),'r') 
hold on 

plot(F1,10*log10(Pxx1 ),'b') 
hold on 

plot(F2,10*log10(Pxx2),'m') 

ylabel('dB re m/s/$$\sqrt{Hz}$$','Interpreter','Latex','Fontsize', 14) 
xlabel('Frequency (Hz)','Fontsize',12) 
title({'Power Spectral Density Estimate'}, 'Fontsize',1 6) 
grid on 

legend('Channel Z','Channel 1','Channel 2') 
figure 

plot(F,10*log10(Pxx),'r') 
hold on 

plot(F1,10*log10(Pxx1),'b') 
hold on 

plot(F2,10*log10(Pxx2),'m') 

ylabel('dB re m/s/$$\sqrt{Hz}$$','Interpreter','Latex','Fontsize', 14) 
xlabel('Frequency (Hz)','Fontsize', 12) 
xlim([0 2]) 

title({'Power Spectral Density Estimate','0 to 2 Hz'}, 'Fontsize', 16) 
grid on 

legend('Channel Z','Channel 1','Channel 2') 

index=abs(cal_data)>=max_amp; 

above_max_amp_time=time(index); 

percentage_above_max_amp_time=(length(above_max_amp_time)/length(time))*100 % will be 
zero with this function 

figure 

histogram(cal_data,400) 

set(gca,'Yscale','log') 
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APPENDIX B. DATA TABLES 


Table 3. YM 39 and YM 38 Standard Deviation and Percentage of 
Time Spent Above Maximum Amplitude Threshold. 


YM 39 

YM38 

Date 

Standard 
deviation (pm/s) 

RMS acceleration 
(pm/s^) 

Percentage above 
maxamplitude 

Date 

Standard 
deviation (pm/s) 

RMS acceleration 
(pm/s') 

Percentage above 
maxamplitude 

13-Jun-08 

191 

155 

17 

13-Jun-08 

0.23 

2.7 

0.00 

14-Jun-08 

143 

65.8 

10 

14-Jun-08 

0.40 

3.6 

0.00 

l-Jul-09 

112 

57.0 

36 

l-Jul-09 

0.45 

3.6 

0.00 

lB-Jul-08 

604 

919 

83 

15-Jul-08 

1.9 

9.5 

0.03 

l-Aug-08 

395 

355 

86 

l-Aug-08 

7.3 

8.1 

0.16 

15-Aug-08 

324 

403 

83 

15-Aug-08 

1.2 

3.7 

0.00 

l-Sep-08 

662 

574 

91 

l-Sep-08 

0.57 

3.8 

0.00 

15-Sep-08 

350 

392 

86 

15-Sep-08 

0.70 

6.1 

0.00 

l-Oct-08 

365 

364 

83 

l-Oct-08 

0.56 

3.9 

0.00 

lB-Oct-08 

375 

303 

77 

15-0ct-08 

0.44 

4.8 

0.00 

l-Nov-08 

399 

358 

80 

l-Nov-08 

1.1 

3.8 

0.01 

15-NOV-08 

219 

279 

74 

15-NOV-08 

1.2 

12.7 

0.00 

13-Dec-08 

200 

489 

42 

13-Dec-08 

1.3 

4.1 

0.01 

14-Dec-08 

277 

303 

33 

14-Dec-08 

1.9 

4.3 

0.04 

l-Jan-09 

727 

185 

84 

l-Jan-09 

O.BB 

4.6 

0.00 

lB-Jan-09 

366 

323 

68 

15-Jan-09 

1.3 

3.9 

0.02 

l-Feb-09 

46.4 

128 

10 

l-Feb-09 

0.42 

3.4 

0.00 

15-Feb-09 

150 

33.0 

30 

lB-Feb-09 

0.43 

3.9 

0.00 

l-Mar-09 

309 

48.4 

57 

l-Mar-09 

1.3 

5.1 

0.00 

lS-Mar-09 

119 

74.6 

54 

15-Mar-09 

0.57 

4.7 

0.00 

l-Apr-09 

96.9 

19.0 

34 

l-Apr-09 

2.9 

13.3 

0.03 

15-Apr-09 

88.0 

71.9 

8 

lB-Apr-09 

1.3 

8.8 

0.01 

l-May-09 

63.9 

86.3 

25 

l-May-09 

0.56 

4.7 

0.00 

15-May-09 

16.1 

28.9 

1 

15-May-09 

0.31 

3.4 

0.00 
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Table 4. 7D J65A and 7D J73A Standard Deviation and Percentage of 
Time Spent Above Maximum Amplitude Threshold. 


7DJ65A 

7DJ73A 

Date 

Standard 
deviation (pm/s) 

RMS acceleration 
(pm/s') 

Percentage above 
max amplitude 

Date 

Standard 
deviation (pm/s) 

RMS acceleration 
(pm/s') 

Percentage above 
max amplitude 

1-Nov-ll 

30.5 

51.7 

9.8 

1-Nov-ll 

28.7 

33.0 

8 

15-Nov-ll 

26.8 

38.5 

6.6 

15-Nov-ll 

24.4 

25.2 

5 

1-Dec-ll 

117 

378 

1.6 

1-Dec-ll 

17.2 

33.9 

1 

15-Dec-ll 

15.3 

38.8 

0.2 

15-Dec-ll 

16.4 

44.0 

0 

l-Jan-12 

29.2 

76.0 

8.8 

l-Jan-12 

166 

391 

2 

15-Jan-12 

31.2 

50.3 

11 

15-Jan-12 

31.5 

40.0 

11 

l-Feb-12 

42.4 

56.9 

24 

l-Feb-12 

39.9 

46.0 

21 

15-Feb-12 

30.5 

48.7 

10 

15-Feb-12 

29.7 

39.4 

9 

l-Mar-12 

40.4 

61.9 

20.0 

l-Mar-12 

34.5 

50.6 

14 

15-Mar-12 

33.3 

73.7 

13.0 

15-Mar-12 

33.4 

74.6 

13 

l-Apr-12 

32.2 

59.3 

12 

l-Apr-12 

38.2 

78.2 

16 

15-Apr-12 

10.0 

25.2 

0.00 

15-Apr-12 

9.2 

21.1 

0 

l-May-12 

25.7 

57.7 

5.2 

l-May-12 

20.8 

51.9 

2 

15-May-12 

14.7 

45.3 

0.1 

15-May-12 

12.2 

31.7 

0 

l-Jun-12 

18.1 

29.4 

1.1 

l-Jun-12 

22.3 

41.3 

3 

15-Jun-12 

11.0 

18.6 

0.00 

15-Jun-12 

10.7 

25.4 

0 

l-Jul-12 

7.2 

29.5 

0.00 

l-Jul-12 

7.6 

31.5 

0 

15-Jul-12 

10 

37 

2.6 

15-Jul-12 

12.5 

35.6 

0 
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Table 5. YM 39 and YM 38 Decibel Values Used in Figure 51. . 


YM 39 

YM 38 

Date 

dBatlHz 

dBat5 Hz 

Date 

dBatlHz 

dBat5 Hz 

13-Jun-08 

-135 

-154 

13-Jun-08 

-147 

-154 

14-Jun-08 


-158 

14-Jun-08 


-153 

l-Jul-09 



l-Jul-09 


-152 

15-Jul-08 

-86 

-103 

15-Jul-08 

-143 

-153 

l-Aug-08 

-97 

-113 

l-Aug-08 

-144 

-148 

15-Aug-08 

-96 

-110 

15-Aug-08 

-144 

-154 

l-Sep-08 

-93 

-108 

l-Sep-08 


-153 

15-Sep-08 

-95 


15-Sep-08 


-149 

l-Oct-08 

-93 

-111 

l-Oct-08 

-141 

-151 

15-Oct-08 

-100 

-116 

15-Oct-08 


-152 

l-Nov-08 

-96 

-111 

l-Nov-08 


-151 

15-NOV-08 

-99 

-113 

15-NOV-08 


-149 

13-Dec-08 

-122 

-137 

13-Dec-08 


-152 

14-Dec-08 

-122 

-126 

14-Dec-08 


-150 

l-Jan-09 



l-Jan-09 


-150 

15-Jan-09 


-114 

15-Jan-09 


-151 

l-Feb-09 

-128 

-151 

l-Feb-09 

-146 

-154 

15-Feb-09 


-143 

15-Feb-09 


-152 

l-Mar-09 



l-Mar-09 


-151 

15-Mar-09 

-122 

-141 

15-Mar-09 

-139 

-148 

l-Apr-09 

-123 

-143 

l-Apr-09 

-143 

-150 

15-Apr-09 

-130 

-132 

15-Apr-09 

-140 

-148 

l-May-09 



l-May-09 


-152 

15-May-09 


-148 

15-May-09 


-153 
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Table 6. 7D J65A and 7D J73A Used in Figure 52. 


7DJ65A 

7DJ73A 

Date 

dBatlHz 

dBat5 Hz 

Date 

dBatlHz 

dBat5 Hz 

1-Nov-ll 

-116 

-151 

1-Nov-ll 

-116 

-155 

15-Nov-ll 

-118 

-154 


-118 

-156 

1-Dec-ll 

-117 

-147 


-118 

-155 

15-Dec-ll 

-112 

-154 

15-Dec-ll 


-155 

l-Jan-12 

-114 

-152 

l-Jan-12 


-153 

15-Jan-12 

-115 

-152 

15-Jan-12 


-149 

l-Feb-12 

-115 

-148 

l-Feb-12 

-114 

-149 

15-Feb-12 

-118 

-150 

15-Feb-12 

-111 

-150 

l-Mar-12 

-114 

-153 

l-Mar-12 

-112 

-146 

15-Mar-12 

-118 

-150 


-109 

-149 

l-Apr-12 

-112 

-147 


-109 

-150 

15-Apr-12 

-121 

-148 

15-Apr-12 

-115 

-156 

l-May-12 

-114 

-153 

l-May-12 

-110 

-150 

15-May-12 

-118 

-150 

15-May-12 

-114 

-153 

l-Jun-12 

-124 

-153 


-117 

-148 

15-Jun-12 

-128 

-152 


-123 

-150 

l-Jul-12 

-121 

-149 

l-Jul-12 

-116 

-149 

15-Jul-12 

-119 

-148 

15-Jul-12 

-108 

-152 
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